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Abstract 

We present two exact implementations of efficient output-sensitive algorithms that compute 
Minkowski sums of two convex polyhedra in M^. We do not assume general position. Namely, 
we handle degenerate input, and produce exact results. We provide a tight bound on the 
exact maximum complexity of Minkowski sums of polytopes in in terms of the number of 
facets of the summand polytopes. The complexity of Minkowski sum structures is directly 
related to the time consumption of our Minkowski sum constructions, as they are output 
sensitive. We demonstrate the effectiveness of our Minkowski-sum constructions through 
simple applications that exploit these operations to detect collision between, and answer 
proximity queries about, two convex polyhedra in M^. 

The algorithms employ variants of a data structure that represents arrangements embed- 
ded on two-dimensional parametric surfaces in M^, and they make use of many operations 
applied to arrangements in these representations. We have developed software components 
that support the arrangement data-structure variants and the operations applied to them. 
These software components are generic, as they can be instantiated with any number type. 
However, our algorithms require only (exact) rational arithmetic. These software components 
together with exact rational-arithmetic enable a robust, efficient, and elegant implementation 
of the Minkowski-sum constructions and the related applications. These software compo- 
nents are provided through a package of the Computational Geometry Algorithm Library 
(Cgal) [5] called Arrangement_on_surf ace_2 [ WFZHOTaj . The code of Cgal in general, 
the Arrangement_on_surf ace_2 package in particular, and all the rest of the code developed 
as part of this thesis adhere to the Generic Programming paradigm and follow the Exact 
Geometric Gomputation paradigm. 

We also present exact implementations of other applications that exploit arrangements 
of arcs of great circles, also known as geodesic arcs, embedded on the sphere. For example, 
we implemented robust polyhedra central-projection and Boolean set-operations applied to 
point sets embedded on the sphere bounded by geodesic arcs. We use them as basic blocks in 
an exact implementation of an efficient algorithm that partitions an assembly of polyhedra 
in M.^ with two hands using infinite translations. This application makes extensive use of 
Minkowski-sum constructions and other operations on arrangements of geodesic arcs embed- 
ded on the sphere. It distinctly shows the importance of exact computation, as imprecise 
computation might result with dismissal of valid partitioning-motions. 

We have produced three movies that explain some of the concepts portrayed in this 
thesis [20] The first movie |FFHL02] explains what Minkowski sums are and demonstrates 
how they are used in various applications. The second movie |FH05j demonstrates the ffist 
method we have developed to construct Minkowski-sums of convex polyhedra. The third 
movie [FSHOSb] illustrates exact construction and maintenance of arrangements induced by 
geodesic arcs and applications that exploit such arrangements. 

Additional information is available at the following web sites: 
http://acg.cs.tau.ac . il/projects/ internal-projects/gaussian-map- cubical 



^Throughout the thesis a number in brackets (e.g., ^Q) refers to the hnk hst starting on page 118, and 
an alphanumeric string in brackets (e.g., [FFHL02] ) is a standard bibhographic reference. 
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|http : //acg. cs . tau. ac . il/projects/ internal-projects/ arrangements-on- surf aces [ 
http : //acg. cs . tau. ac . il/projects/internal-projects/exact-complexity-of - 
minkowski-sums, and 

[http : //acg. cs . tau. ac . il/projects/ internal-projects/ arr-geodesic-spherel 

Auxiliary programs, source code, data sets, and documentation can be downloaded from 
http : //www. cs . tau. ac . il/~ef if /Minkowski-sum. 
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Introduction 



Let P and Q be two closed convex polyhedra in W^. The Minkowski sum of P and Q is the 
convex polyhedron M = P (B Q = {p + q \ p & P, q & Q} ■ A polyhedron P translated by a 
vector t is denoted by P*. Collision Detection is a procedure that determines whether P and 
Q overlap. The Separation Distance 7i{P, Q) and the Penetration Depth 6{P, Q) defined as 

7r{P,Q) =min{||t|| |P*nQ ^ 0,t G M'^} , 
6{P,Q) = inf{||t|| |P*ng = 0,t G M'^} 

are the minimum distances by which P has to be translated so that P and Q intersect or 
become interior disjoint, respectively. The problems of finding the distances above can also 
be posed given a normalized vector r that represents a direction, in which case the minimum 
distance sought is in direction r. The Directional Penetration-Depth, for example, is defined 

as 

6r{P, Q) = inf {a | P"^^ n Q = 0, a G M} . 

1.1 Main Contribution 

We present two exact and robust implementations of efficient output-sensitive algorithms 
to compute the Minkowski sum |FFHL02] of two convex polyhedra, polytopes for short, in 
|FH05l [FH071 IFSHOSal lBFH+09aj . The implementations are exact and robust, as they 
handle all degenerate cases, and guarantee exact results. We demonstrate the effectiveness 
of our Minkowski-sum computations through simple applications that exploit these opera- 
tions to detect collision, and compute the Euclidean separation-distance between, and the 
directional penetration depth of, two polytopes in R^. 
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We compare our Minkowski-sum constructions with three other methods that produce 
exact resuhs we are aware of. One is a simple method that computes the convex hull of the 
pairwise sums of vertices of two polytopes. The second is based on Nef polyhedra embedded 
on the sphere, and the third is based on linear programming. We conducted experiments 
with a broad family of polytopes and compared the performance of our methods with the 
performance of the other. The results reported in Table 13.51 clearly shows that both our 
methods are significantly faster. 

Each method we have developed uses a different variant of Gaussian maps, also known 
as normal diagrams or slope diagrams, to maintain dual representations of polytopes. Each 
method employs a different variant of a generic data-structure that represents arrangement 
embedded on two-dimensional parametric surfaces in to maintain the dual representations. 
(Arrangements embedded on two-dimensional parametric surfaces are subdivisions of the 
surface, as induced by curves embedded on the surface. They play a central role in this 
thesis; see Chapter |2] for a formal definition and Figure [?!T] for an illustration.) Each method 
makes use of many operations applied to arrangements in the corresponding representations. 

We present a tight bound on the exact maximum complexity of Minkowski sums of 
polytopes in M'^ |FHW07j . In particular, we prove that the maximum number of facets 
of the Minkowski sum of k polytopes with mi,m2, . . . ,mk facets respectively is bounded 
from above by X]i<j<i<fc(^"^* ~ 5){2mj — 5) + X]i<j<A;"^j + (2)- Given k positive integers 
mi, m2, . . . , mfc, we describe how to construct k polytopes with corresponding number of 
facets, such that the number of facets of their Minkowski sum is exactly X]i<i<j<fc(2^i " 
5)(2mj — 5) + X]i<j<A: + (2)- When k = 2, for example, the expression above reduces to 
4mim2 — 9mi — 9m2 + 26. 

We also present an exact implementation of an efficient algorithm that partitions an as- 
sembly of polyhedra in with two hands using infinite translations |FH08j . This application 
makes extensive use of Minkowski-sum constructions and other operations on arrangements 
of arcs of great circles, also known as geodesic arcs, embedded on the sphere, such as poly- 
hedra central-projection and Boolean set-operations of point sets embedded on the sphere 
bounded by geodesic arcs. The assembly partitioning demonstrates the importance of exact 
computation, as imprecise computation might result with dismissal of valid partitioning- 
motions. 

Both methods that construct Minkowski sums and the related applications are imple- 
mented on top of the Computational Geometry Algorithms Library (Cgal) )FT07j . and 
are mainly based on the arrang ement package of the library |FWH04l IWFZHOTbi lFHK+07[ 
lBFH+07l[B"FH+09bj . which supports arrangements embedded on two-dimensional paramet- 
ric surfaces and operations on them. We have redesigned, re-implemented, and significantly 
extended the package exploiting advanced programming techniques to yield a package that 
is easy to use, to extend, and to adapt to a variety of applications. 

The package contains a general framework for computing the zone of a single curve em- 
bedded on a two-dimensional parametric surface and a general framework for sweeping a set 
of such curves. The former framework is used, for example, to insert curves one at a time into 
the arrangement. The latter framework is used, for example, to compute the arrangement 
induced by a collection of curves, and to compute the overlay of two arrangements. Other 
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operations, such as point location and vertical ray shooting, are supported as well. 

The design dictates the separation between the topological and geometric aspects of 
the two-dimensional subdivision. This separation is advantageous, as it allows users to 
employ the package with their own representation of any special family of curves. They 
must however supply a relevant component called geometry traits class that handles the 
specific family of curves they are interested in, which mainly involves algebraic computation. 
The separation is enabled by a modular design and conveniently implemented within the 
generic-programming paradigm. The separation is a key aspect of the package, as well as of 
other central Cgal components, such as the various triangulation packages [BDTYOO] and 
convex- hull algorithms (see |cga07| for more details), has been forced since its early stages, 
and heightened by our new design. 

The package comes with many geometric traits classes that handle all kinds of curves 
organized in a structured hierarchy. In particular, we mention the geometric traits class that 
handles continuous piecewise linear curves, referred to as polylines and the geometric traits 
class that handles geodesic arcs embedded on the sphere [FSHnSbl IFSHOSai lBFH+n9aj . The 
former are of particular interest, as they can be used to approximate more complex curves 
in the plane. At the same time they are easier to deal with in comparison to higher-degree 
algebraic curves, as rational arithmetic is sufficient to carry out exact computations on 
polylines. The latter is broadly used by the assembly-partitioning application and by the 
second method that constructs Minkowski sums. 

The implementation of the package as well as the implementation our Minkowski-sum 
constructions, collision detection, and assembly partitioning applications handle degener- 
ate input and produce exact results, as long as the underlying number type supports the 
arithmetic operations -|-, — , *, and / in unlimited precision over the rationals,^ such as the 
rational number type CGAL: :Gmpq based on GMP — Gnu's Multi Precision library [T^ . 

1.2 Background: Minkowski Sums 

Minkowski sums are closely related to proximity queries |LM04] . For example, the mini- 
mum separation distance between two polytopes P and Q is equal to the minimum distance 
between the origin and the boundary of the Minkowski sum of P and the reflection of Q 
through the origin |CC86j . Computing Minkowski sums, collision detection and proximity 
calculation constitute fundamental tasks in computational geometry [HKL04|, ILM041 ISha04j . 
These operations are ubiquitous in robotics, solid modeling, design automation, manufac- 
turing, assembly planning, virtual prototyping, and many more domains; see, e.g., |BdLT97[ 
IEUR921 [KR9T1 iLatOTl IVKSMOSj . 

The wide spectrum of ideas expressed in the massive amount of literature published about 
the subject during the last three decades has inspired the development of quite a few useful 
solutions. For an extensive overview about the subject and a comprehensive list of packages 
see |LM04] . However, only recent advances in the implementation of computational-geometry 
algorithms and data structures made our exact, robust, and efficient implementation possi- 

^ Commonly referred to as a field number type. 
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ble. The exact geometric-computation paradigm |Yap04| designed for implementing compu- 
tational geometry algorithms particularly prevails in Cgal. While in general the underlying 
arithmetic is highly time consuming compared to machine floating-point arithmetic, major 
efficiency is gained by computing predicates only to sufficient precision to evaluate them 
correctly; see Section [1. 3. 21 and e.g., |PF06] . 

Various methods to compute the Minkowski sum of two polyhedra in have been 
proposed. The goal is typically to compute the boundary of the sum provided in some 
representation. The combinatorial complexity of the Minkowski sum of two polyhedra of m 
and n features respectively can be as high as G(m^n^). One common approach to compute 
it is to decompose each polyhedron into convex pieces, compute pairwise Minkowski sums of 
pieces of the two, and finally the union of the pairwise sums. Computing the exact Minkowski 
sum of non-convex polyhedra is naturally expensive. Therefore, researchers have focused on 
computing an approximation that satisfies some criteria, such as the algorithm presented by 
Varadhan and Manocha [VM06j . They guarantee a two-sides Hausdorff distance bound on 
the approximation, and ensure that it has the same number of connected components as the 
exact Minkowski sum. Computing the Minkowski sum of two convex polyhedra remains a 
key operation, and this is what we focus on. The combinatorial complexity of the sum can be 
as high as 0(mn) when both polyhedra are convex. For the complexity of the intermediate 
case, where only one polyhedron is convex, cf. |AS97[ [S"ha04j . 

Convex decomposition is not always possible, as in the presence of non-convex curved 
objects. In these cases other techniques must be applied, such as approximations using 
polynomial/rational curves in 2D [kLsKE98] . Seong at al. |SKS02] proposed an algorithm 
to compute Minkowski sums of a subclass of objects; that is, surfaces generated by slope- 
monotone closed curves. Plato and Halperin [AFH02j presented algorithms based on Cgal 
for robust construction of planar Minkowski sums. While the citations in this paragraph 
refer to computations of Minkowski sums of non-convex polyhedra, and we concentrate on 
the convex cases, the latter is of particular interest, as our method makes heavy use of 
the same software components, in particular the Cgal arrangement package, which went 
through a few phases of improvements |FWH04l IWFZH07bl lBFH+07[ lBFH+09bj since its 
usage in [AFH02j : see Section [1.3.31 for more details. 

A particular accomplishment of the kinetic framework in two dimensions introduced by 
Guibas et al. |GRS83j was the definition of the convolution operation in two dimensions, a 
superset of the Minkowski sum operation, and its exploitation in a variety of algorithmic 
problems. Basch et al. extended its predecessor concepts and presented an algorithm to 
compute the convolution in three dimensions |BGRR96] . An output-sensitive algorithm 
for computing Minkowski sums of polytopes was introduced in |GS87j . Gritzmann and 
Sturmfels [GS93j obtained a polynomial time algorithm in the input and output sizes for 
computing Minkowski sums of k polytopes in for a fixed dimension rf, and Fukuda |Fuk04] 
provided an output-sensitive polynomial algorithm for variable d and k. Ghosh |Gho93] 
presented a unified algorithm for computing 2D and 3D Minkowski sums of both convex and 
non-convex polyhedra based on a slope diagram representation. Computing the Minkowski 
sum amounts to computing the slope diagrams of the two objects, merging them (see details 
in Section I3.2.2[ ) and extracting the boundary of the Minkowski sum from the merged 
diagram. Wu et al. |WSD03j introduced an improved version of Ghosh' algorithm for convex 
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polyhedra using vector operations. Bekker and Roerdink [BROl] provided a variation on the 
same idea. The slope diagram of a 3D convex polyhedron can be represented as a 2D object, 
essentially reducing the problem to a lower dimension. We follow the same approach, but 
use exact computation. 

We postpone a formal definition of arrangements to the next chapter. In fact, we tem- 
porary put Minkowski sums and arrangements aside to provide relevant programming back- 
ground material. 



1.3 Background: Programming 
1.3.1 Generic Programming 

Several definitions of the term generic programming have been proposed since it was first 
coined around the early sixties, along with the introduction of the Lisp programming lan- 
guage. Here we confine ourself to the classic notion first described by Musser et al. |MS88j . 
who considered generic programming as a discipline that consists of the gradual lifting of 
concrete algorithms abstracting over details, while retaining the algorithm semantics and 
efficiency. Within this context, several approaches have been put into trial through the in- 
troduction of new features in existing computer languages, or even new computer languages 
all together. The software described in this thesis is written in C++, a programming language 
that is well equipped for writing software according to the generic-programming paradigm 
through the extensive use of class templates and function templates. 

Concepts and Models 

One crucial abstraction supported by all contemporary computer languages is the subroutine 
(also known as procedure or function, depending on the programming language). Another 
abstraction supported by C++ is that of abstract data typing, where a new data type is 
defined together with its basic operations. C++ also supports object-oriented programming, 
which emphasizes packaging data and functionality together into units within a running 
program, and is manifested in hierarchies of polymorphic data-types related by inheritance. 
It allows referring to a value and manipulating it without needing to specify its exact type. As 
a consequence, one can write a single function that operates on a number of types within an 
inheritance hierarchy. Generic programming identifies a more powerful abstraction (perhaps 
less tangible than other abstractions). It is a formal hierarchy of polymorphic abstract 
requirements on data types referred to as concepts, and a set of classes that conform precisely 
to the specified requirements, referred to as models. Models that describe behaviors are 
referred to as traits classes |Mye98| . Traits classes typically add a level of indirection in 
template instantiation to avoid accreting parameters to templates. 

A generic algorithm has two parts: The actual instructions that describe the steps of the 
algorithm, and a set of requirements that specify which properties its argument types must 
satisfy. The following swap ( ) function is an example of the first part of a generic algorithm. 
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template <typenaine T> void swap(T & a, T & b) { 
T tmp = a; a = b; b = tmp; 

_} 

When the function call is compiled, it is instantiated with a data type that must have an 
assignment operator. A data type that fulfils this requirement is a model of a concept 
commonly called Assignable |Aus99j . The int data type, for example, is a model of this 
concept, so it can be used to instantiate the function template |Aus99j [23] . 

A concept is a set of requirements divided into four categories, namely, associated types, 
valid expressions, invariants, and complexity guarantees. When a type meets all requirements 
of a concept, the type is considered a model of the concept. When a concept extends the 
requirements of another concept, the former is said to be a refinement of the latter. 

Associated Types are auxiliary types. For example, a type that represents a two-dimensional 
point, namely Point_2, is required by the arrangement geometry traits concept; see 
Section 12.51 

Valid Expressions are C++ expressions that must compile successfully. For example, p = 
q, where p and q are both objects of type Point_2. Valid expressions identify the set 
of operations a model of the concept must be able to perform. 

Invariants are run-time characteristics such as time and space complexity bounds. In our 
context invariants typically take the form of preconditions and postconditions, which 
must always be satisfied. For example, a condition that requires that an input point 
p lies on an input curve c on invocation to a predicate that accepts both p and c as 
parameters. Having preconditions typically minimizes the concept, as the operations 
provided by a model must operate only on restricted arguments. Formally, removing 
preconditions from, and introducing postconditions to, a requirement set results with 
a refined concept. 

Complexity Guarantees are maximum limits on the computing resources consumed by 
the various expressions. 



Traits Classes 



The name "traits class" comes from a standard C++ design pattern |Mye98| , which provides 
a way of associating information with a compile-time entity (typically a type). For example, 
the standard class-template std: : iterator_traits<T> looks roughly like this: 
template <typename Iterator> struct iterator_traits { 



typedef 
typedef 
typedef 
typedef 
typedef 



iterator_category ; 

value_type ; 

dif f erence_type ; 

pointer; 

reference ; 



Iterators play an important role in generic programming: A function that operates on a 
range of objects usually accepts two iterators that specify this range. The traits' value_type 
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specifies the type of object the iterators are pointing at, while the iterator_category can 
be used to select more efficient algorithms depending on the iterator's capabilities. 

A key property of trait classes is that they are non-intrusive. Namely, they allow us to 
associate information with arbitrary types, without interfering with the internal representa- 
tion of those types. Thus, it is possible to define a traits class also for built-in types and 
types defined in third-party libraries. 

Within the context of Cgal, for example, a typical traits class is required to define nested 
types of geometric objects and support predicates involving objects of these types. Some 
algorithms also require the provision of constructions by the traits class. 

Let us continue with an easy geometric example. Consider a function that accepts a 
set of points, given by the range [pts_begin, pts_end),^ and computes the minimal axis- 
parallel rectangle that contains all points in the range. It does so by locating the points with 
extremal x and y-coordinates, and then constructs the bounding iso-rectangle accordingly: 

template <typename Input Iterator , typename Traits> 
typename Traits :: Iso_rectangle_2 

bounding_rectangle(InputIterator pts_begin, Inputlterator pts_end) { 
Traits traits; 

Inputlterator curr = pts_begin; 
Inputlterator left, right, top, bottom; 
left = right = top = bottom = curr++; 
while (curr++ != pts_end) { 

if (traits . compare_x(*curr, *left) == SMALLER) left = curr; 

else if (traits . compare_x(*curr, *right) == LARGER) right = curr; 

if (traits . compare_y(*curr, *bottom) == SMALLER) bottom = curr; 

else if (traits . compare_y(*curr, *top) == LARGER) top = curr; 

} 

return traits. construct_iso_rectangle(*left, *right, *bottom, *top) ; 

_A 

The requirements that an instantiated traits class must satisfy in this case are as fol- 
lows: It has to defined the nested type Iso_rectangle_2 (and implicitly also a point type 
say Point_2). Moreover, it should supply two three- valued predicates^ that compare two 
points by their ^-coordinates and by their ^/-coordinates, respectively. It should also support 
the construction of an axis-parallel iso-rectangle from four points that specify its extremal 
X and y-coordinates. Note, however, that the actual representation of points and rectan- 
gles (the coordinate system, the number-type used to represent the coordinates, etc.) and 
the implementation of the traits-class operations is entirely decoupled from the function 
bounding_rectcingle() we have introduced. 

Consider an imaginary generic implementation of a data structure that handles ge- 
ometric arrangements embedded on two-dimensional parametric surfaces in space called 
Arrangement_on_surf ace_2. Its prototype is listed below. This template class must be 



^This notation means that pts_begin points to the first point in the range, while pts_end points after 
the range ends (it is therefore called a past-the-end iterator, and need not point to any valid point object). 
^The predicate return value is SMALLER, EQUAL, or LARGER. 
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instantiated with a traits class that in turn defines a type that represents a certain family of 
curves, and some functions (or function objects; see |T0] for an exact definition) that operate 
on curves of this family. 

template <typename Traits> class Arrangement_on_surf ace_2 { ... } ; 

Traits classes that handle families of curves embedded on parametric surfaces are intricate, as 
they model involved concepts. The precise definitions of these concepts and their refinement 
hierarchy are described in Section 12.51 

One important objective is to minimize the set of requirements the traits concept imposes. 
A tight traits concept may save tremendously in analysis and programming of classes that 
model the concept. Another important reason for striving for the minimal set of requirements 
is to avoid computing the same algebraic entity in different ways. The importance of this 
is amplified in the context of computational geometry, as a non tight model that consists of 
duplicate, but slightly different, implementations of the same algebraic entity, can lead to 
artificial degenerate conditions, which in turn can drastically impair the performance. 

Most traits classes in Cgal are parameterized by a model of the Kernel concept. The 
choice of a particular model determines, among the other, the type of arithmetic used, as 
explained in the following sections. One can easily switch between different models of the 
Kernel concept, but here lies a trap, as Section [1.3.21 reveals. A kernel model that supports 
exact arithmetic must be used to ensure robustness, although inexact arithmetic could be 
used at a certain risk. 

Libraries 

Alexander Stepanov began exploring the potential of compile-time polymorphism for rev- 
olutionizing software development in 1979. With the help of several other researchers his 
work evolved into a prime generic-programming library — the Standard Template Library 
(STL). This library had became part of the C++ standard library in 1994, approximately 
one year before early development of Cgal started; see Section 11.3.31 for details about the 
evolution of Cgal. 

Through the years a few other generic-programming libraries emerged. One notable li- 
brary in the context of computational geometry is Leda (Library of Efficient Data Types and 
Algorithms), a library of combinatorial and geometric data types and algorithms jMNOO] |17j . 
Early development of Leda started in 1988, ten years before the first public release of Cgal 
became available. In some sense Leda is a predecessor of Cgal, although the two libraries 
are headed in different directions. While Leda is mostly a large collection of fundamental 
graph related and general purpose data-structures and algorithms, Cgal is a collection of 
large and complex data-structures and algorithms focusing on geometry. 

A noticeable influence on generic programming is conducted by the Boost online commu- 
nity, which encourages the development of free C++ software gathered in the Boost library 
collection pj. It is a large set of portable and high quality C++ libraries that work well 
with, and are in the same spirit as, the C++ Standard Template Library. The Boost Graph 
Library (BGL) |SLL02j . which consists of generic graph-algorithms, serves a particularly 
important role in our context. An arrangement instance, for example, can be adapted as a 
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BGL graph, and passed as input to generic algorithms aheady implemented in the BGL, 
such as the Dijkstra shortest path algorithm. We use the BGL to compute the strongly con- 
nected components of a directed graph — a phase in the assembly partitioning operation; 
see Section 15.3.81 for more details about the application of this operation. 



1.3.2 Geometric Programming 

Implementing geometric algorithms and data structures is notoriously difficult, much harder 
than may seem when just considering the algorithm as described in a paper or a book. 
In the traditional computational-geometry literature two assumptions are usually made to 
simplify the design and analysis of geometric algorithms. First, inputs are in "general po- 
sition". That is, degenerate input (e.g., three curves intersecting at a common point) is 
precluded. Secondly, operations on real numbers yield accurate results (the "real Ram" 
model [PS85j . which also assumes that each basic operation takes constant time). Unfortu- 
nately, these assumptions do not hold in practice, as degenerate input is commonplace in 
practical applications and numerical errors are inevitable. Thus, an algorithm implemented 
without keeping this in mind may yield incorrect results, get into an infinite loop, or just 
crash, while running on a degenerate, or nearly degenerate, input (see |KMP"'"08i ISchOO] for 
examples). These pitfalls have become well known, and have been the subject of intensive 
research |SchOO| |Yap04| . 



Indeed, the last decade has seen significant progress in the development of software for 
computational geometry. The mission of such a task, which Kettner and Naher |KN04] 
call geometric programming, is to develop software that is correct, efficient, fiexible (namely 
adaptable and extensible'*), and easy to use. 



Separation of Topology and Geometry 

The use of the generic-programming paradigm enables a convenient separation between the 
topology and the geometry of data structures.^ This is a key aspect in the design of ge- 
ometric software, and is put into practice, for example, in the design of Cgal polyhedra, 
Cgal triangulations, and our Cgal arrangements. This separation allows the convenient 
abstraction of algorithms and data structures in combinatorial and topological terms, re- 
gardless of the specific geometry of the objects at hand and the algebra used to represent 
them. This abstraction is realized through class and function templates that represent spe- 
cific data-structures and algorithmic frameworks, respectively. Consider again our imaginary 
Arrangement_on_surf ace_2 class template from the previous section; its improved prototype 
is listed below. It is instantiated with two classes. The first, referred to as a geometric traits 
class, defines the set of geometric-object types and the operations on objects of these types. 
The second, defines the topological-object types and the operations required to maintain the 

'^Adaptability refers to the ability to incorporate existing user code, and extendibility refers to the abihty 
to enhance the software with more code in the same style. 

^In this context, we sometimes say combinatorics instead of topology, and say algebra or numerics instead 
of geometry. We always mean the same thing: The separation between the abstract, graph-like structure 
(the topology) from the actual embedding on the surface (the geometry). 
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incident relations among objects of these types, 
template <typename Geometry_traits , typename Topology _traits> 
class Arrangement_on_surf ace_2 { ... } ; 

An immediate advantage of the separation between the topology and the geometry of data 
structures is that users with limited expertise in computational geometry can employ the 
data structure with their own special type of objects. They must however supply the relevant 
traits class, which mainly involve algebraic computations. A traits class also encapsulates the 
number types used to represent coordinates of geometric objects and to carry out algebraic 
operations on them. It encapsulates the type of coordinate system used (e.g., Cartesian, 
Homogeneous), and the geometric or algebraic computation methods themselves. Naturally, 
a prospective user of the package that develops a traits class would like to face as few 
requirements as possible in terms of traits development. 

Another advantage gained by the use of generic programming is the convenient handling 
of numerical issues to expedite exact geometric computation. We arrive at this conclusion at 
the end of the next section. In a geometric algorithm each computational step is either a con- 
struction step or a conditional step based on the result of a predicate. The former produces 
a new geometric object such as the intersection point of two segments. The latter typically 
computes the sign of an expression used by the program control. Different computational 
paths lead to results with different combinatorial characteristics. Although numerical errors 
can sometimes be tolerated and interpreted as small perturbations in the input, they may 
lead to invalid combinatorial structures or inconsistent state during a program execution. 
Thus, it suffices to ensure that all predicates are evaluated correctly to eliminate inconsis- 
tencies and guarantee combinatorially correct results. This is easier said than done, but 
nowadays possible, as the next section exposes. 

Exact Geometric Computation 

The need for robust software implementations of computational-geometry algorithms has 
driven many researchers over the last decade to develop variants of the classic algorithms that 
are less susceptible to degenerate inputs. The approaches taken to overcome the difficulties 
in robustly implementing geometric algorithms can be roughly divided into two categories: 
(i) Exact computing, and (ii) fixed-precision approximation. In the latter approach the 
algorithms are modified so that they can consistently cope with the limited precision of 
computer arithmetic. In the former, which is the approach taken by Cgal in general and the 
arrangement package in particular, ideal computer arithmetic is emulated for the specific type 
of objects being manipulated, and the code is prepared for successfully handling degenerate 
input. 

Advances in computer algebra enabled the development of efficient software libraries 
that offer exact arithmetic manipulations on unbounded integers, rational numbers (GMP 
— Gnu's multi-precision library [12]), and algebraic numbers (the Core library [KLPY99] [S] 
and the numerical facilities of Leda [MNOOt Chapter 4] [17]). These exact-number types 
serve as fundamental building-blocks in the robust implementation of many geometric ap- 
plications in general (see |Yap04| for a review) and of those that employ arrangements in 
particular. 
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Exact Geometric Computation (EGC), as summarized by Yap |Yap04| , simply amounts to 
ensuring that we never err in predicate evaluations. EGC represents a significant relaxation 
from the naive concept of numerical exactness. We only need to compute to sufficient 
precision to make the correct predicate evaluation. This has led to the development of several 
techniques such as precision-driven computation, lazy evaluation, adaptive computation, and 
floating-point filters, some of which are implemented in Cgal, such as numerical filtering. 
Here, computation is carried out using a number type that supports only inexact arithmetic 
(e.g., double-precision fioating-point arithmetic), while applying a filter that checks whether 
the computation has reached a stage of uncertainty, an event referred to as a filter failure 
in the hacker's jargon. When a filter failure occurs, the computation is re-done using exact 
arithmetic. 

Switching between number types and exact computation techniques, and choosing the 
appropriate components that best suit the application needs, is conveniently enabled through 
the generic-programming paradigm, as it typically requires only a minor code change refiected 
in the instantiating of just a few data types. 



1.3.3 Computational Geometry Algorithms Library 

Cgal Chronicles 

Several research groups in Europe started to develop small geometry libraries on their own 
in the early 1990s. A consortium of several sites in Europe and Israel was founded in 
1995 to cultivate the labor of these groups and gather their produce in a common library 
called Cgal — the Computational Geometry Algorithms Library. The goal was to promote 
the research in computational geometry and translate the results into useful, reliable, and 
efficient programs for industrial and academic applications |Ove96i |cgaU7[ IKN041 [FGK"'"00j . 
the very same goal that governs Cgal development efforts to date. 

An INRIA startup. Geometry Factory [UJ was founded in January 2003. The com- 
pany sells Cgal commercial licenses, support for Cgal, and customized developments based 
on Cgal. 

In November 2003, when Version 3.0 was released, Cgal became an Open Source 
Project [5], allowing new contributions from various resources. Most parts of Cgal are 
now distributed under the GNU Lesser General Public License (GNU LGPL). 

Cgal has evolved through the years and is now representing the state-of-art in imple- 
menting computational geometry software in many areas. The implementations of the Cgal 
software modules described in this thesis are complete and robust, as they handle all de- 
generate cases. They rigorously adhere to the generic-programming paradigm to overcome 
problems encountered when effective computational geometry software is implemented. A 
glimpse at the structure of Cgal is given in the following subsection. 



Cgal Content 



Cgal is written in C++ according to the generic-programming paradigm described above. 
It has a common programming style, which is very similar to that of the STL. Its ap- 
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plication programming-interface (API) is homogeneous, and allows for a convenient and 
consistent interfacing with other software packages and applications. The library consists 
of about 900,000 lines of code divided among approximately 4,000 files. Cgal also comes 
with numerous examples and demos. The manual comprises about 3,500 pages. There are 
approximately 65 chapters arranged in 14 parts. The Arrangement_on_surf ace_2 package, 
for example, consists of about 140,000 lines of code divided among approximately 300 files, 
described in about 300 pages of a didactic manual. 

One distinguished piece consists of the geometric kernels |FGK+00j . A geometric kernel 
consist of types of constant size non- modifiable geometric primitive objects (e.g., points, 
lines, triangles, circles, etc.) and operations on objects of these types. 

Another distinguished piece, referred to as the "Support Library" consists of non-geometric 
facilities, such as circulators, random generators, and I/O support for interfacing Cgal with 
various visualization tools (i.e., input and output streams). An important contribution of 
this piece is the number- type support. This piece also contains extensive debugging utilities 
that handle warnings and errors that may result from unfulfilled conditions. 

The rest of the library offers a collection of geometric data structures and algorithms such 
as convex hull, polygons and polyhedra and operations on them (Boolean operations, polygon 
offsetting), 2D and 3D triangulations, Voronoi diagrams, surface meshing and surface subdi- 
vision, search structures, geometric optimization, interpolation, and kinetic data-structures. 
The 2D arrangements and its related data-structures naturally fit in. These data structures 
and algorithms are parameterized by traits classes that define the interface between them 
and the primitives they use. In many cases, a kernel can be used as a traits class, or at least 
the subtypes of a kernel can be used as components of traits classes for these data structures 
and algorithms. 



Cgal Arrangement Package History 

Cgal contains an elaborate and efficient implementation of a generic data-structure that 
represents an arrangement induced by general types of curves embedded on a two-dimensional 
parametric surface in M^, but it has not been like this from the beginning. The first version 
of the Cgal arrangement package supported only line segments, circular arcs, and restricted 
types of parabolas embedded in the plane. 

While the first version supported only limited types of curves, it was originally designed 
with the vision of supporting general curves |FHH+00[ IHanOOl [HHOO] . This vision was 
refiected, among the other, through the separation between the topological and the geometric 
aspects enabled by the generic-programming paradigm (see Section ll.3.2p . Most of the 
principles related to the topology, e.g., the use of a doubly- connected edge list (Dcel) to 
maintain the incident relations between the arrangement cells (i.e., vertices, halfedges, and 
faces) were conceived from the start. However, the types of curves that induce arrangements 
(see Section 12.51) gradually expanded. A couple of years after the introduction of Cgal 
arrangement Wein extended its implementation to support arcs of conies |Wei02] . The 
arsenal of geometric traits continues to grow to date. 
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Following the requirements that emerged from the ECG project^, together with Wein 
we improved and refined the software design of the Cgal arrangement package [FWH04| . 
This new design formed a common platform for a preliminary comparison between different 
approaches to handle arcs of conies [FHW"'"04] . Two years later in a joint effort with Wein 
and Zukerman the whole package was revamped |WFZH05] leading to more compact, easier- 
to-use, and efficient code. 

Cgal Version 3.2 released in 2006 included an arrangement package that supported only 
bounded curves in the plane. This forced users to clip unbounded curves before inserting 
them into the arrangement; it was the user responsibility to clip without loss of information. 
However, this solution is generally inconvenient and outright insufficient for some applica- 
tions. For example, representing the minimization diagram defined by the lower envelope of 
unbounded surfaces in |Mey06| generally requires more than one unbounded face, whereas 
an arrangement of bounded clipped curves contains a single unbounded face. 

Cgal Version 3.3 released a year later in 2007 already included an arrangement pack- 
age that handled unbounded planer curves. As a matter of fact it included much more. 
Together with Berberich, Melhorn, and Wein we observed the possibility to maximize code 
reuse by generalizing the various algorithms applied to arrangements and their implementa- 
tions so that they could be employed on a large class of surfaces and curves embedded on 
them |BFH+07[ lBFH"'"09b] . Indeed, the algorithms and their implementations were designed 
with the vision of supporting general curves embedded on parametric surfaces. However, only 
a few geometric traits-classes that supported unbounded curves in the plane were included 
in Version 3.3. 

A future version of Cgal, expected to be released in 2010, is planned to include an ar- 
rangement package that constructs, maintains, and operates on arrangements embedded on 
certain two-dimensional orientable parametric surfaces. The package already exists as a pro- 
totypical Cgal package under the new name Arrangement_on_surf ace_2 to better refiect 
its capabilities. For example, it includes (i) a geometric traits that handles geodesic arcs em- 
bedded on the sphere |FSH08b[ iBFH+OQaj . (ii) a geometric traits that handles intersections 
between quadric surfaces embedded on a quadric |BFH+07l [BFH+09aj . and (iii) a geometric 
traits that handles intersections between arbitrary algebraic surfaces and a parameterized 
Dupin cyclide embedded on the Dupin cyclide |BK08l iBFH+OQaj . The references to the 
arrangement software in this thesis in general and in Chapter [2] in particular pertain to this 
latest version. 

The leap in arrangement technology expressed by the ability to construct and maintain 
arrangements embedded on two-dimensional parametric surfaces immediately affects other 
components in Cgal, such as the Boolean_set_operations_2 package. Only little effort 
is now required to support Boolean set-operations on point sets bounded by general curves 
embedded on two-dimensional parametric surfaces. 



^ECG is a Shared-Cost RTD (FET Open) Project of the European Union devoted to Effective Compu- 
tational Geometry for curves and surfaces [Tj. 
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1.4 Thesis Outline and Related Publications 

The rest of this thesis is organized as follows. In Chapter [2] we give an overview of the Cgal 
arrangement package, which provides the common infrastructure for all software solutions 
described in this thesis. This chapter contains selected sections from several papers and 
from a book chapter we have co-authored, and it provides new material that has not been 
published yet. In particular, the chapter contains a description of the architecture of the ar- 
rangement generic data-structure, part of which also appeared in the chapter Arrangement 
of the book Effective Computational Ceometry for Curves and Surfaces |FHK"'"07] . The 
chapter contains a description of advanced programming techniques applied in the context 
of arrangement, parts of which appeared in a joint work with R. Wein and B. Zukererman, 
and published in the journal Computational Ceometry — Theory and Application (special 
issue on Cgal) |WFZH07b] . Preliminary results were first introduced in (i) the proceed- 
ings of the 12*^^ Annual European Symposium on Algorithms (ESA) |FWH04j and (ii) the 
1^* Workshop of Library- Centric Software Design |WFZH05] . The chapter also presents ar- 
rangements induced by general curves and embedded on parametric surfaces, and it offers 
a detailed description of a particular type of arrangement, namely arrangements of geodesic 
arcs on the sphere. The presentation of general arrangements embedded on parametric sur- 
faces is based on a joint work with E. Berberich, K. Melhorn, and R. Wein. Preliminary 
results of this work were first published at the 23'''^ European Workshop on Computational 
Ceometry (EWCC) |BFHW07] . Improved results appeared in the proceeding of the 15*^ 
Annual European Symposium on Algorithms (ESA) |BFH+07] . and mature results were pre- 
sented in a manuscript [BFH+OQb] . The discussion about particular arrangements embedded 
on spheres is based on joint work with O. Setter. Preliminary results of this work were first 
published at the 24*'' European Workshop on Computational Ceometry (EWCC) |FSH08a] . 
A movie rendering the results of this work was presented at the 24*^^ Annual Symposium on 
Computational Ceometry (SoCC). The proceeding of this symposium contains the related 
extended abstract [FSHOSbj . Mature results were presented in a manuscript [BFH+OQa) . 

In Chapter [3] we present two different complete, yet efficient, implementations of output- 
sensitive algorithms to compute the Minkowski sum of two polytopes in M? . We describe 
how the input polytopes in polyhedral-mesh representation both methods accept are con- 
verted into the corresponding internal representations unique to each method. We provide 
the theoretical concepts both methods rely on, and detailed descriptions specific to each 
method. The first method was also published in the journal Computer Aided Design |FH07j . 
A preliminary version of this paper appeared in the proceedings of the 8*'^ Workshop on 
Algorithm Engineering and Experimentation (Alenex'06) |FH06] . The second method ap- 
peared in the publications |FSH08a| IFSHOSbl iBFH+OQaj mentioned above. We compare our 
Minkowski-sum constructions with other methods that produce exact results, and provide a 
summary of a performance analysis of our methods. 

Chapter H] provides a tight bound on the exact maximum complexity of Minkowski sums 
of k polytopes in in terms of the number of facets of the polytope summands. It is based 
on collaborative work with C. Weibel. The results of this work were introduced in the pro- 
ceedings of the 23'"'^ annual Symposium on Computational Ceometry (SoCC) [FHW07] . and 
were accepted for publication in the journal Discrete and Computational Ceometry [FHWj . 
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We use this opportunity to thank Shakhar Smorodinsky for fruitful discussions conducted 
while we were investigating the bound above. The chapter provides a proof of the upper 
bound, and establishes the lower bound through a construction procedure. 

Chapter O introduces an exact implementation of an efficient algorithm to obtain a parti- 
tioning motion given an assembly of polyhedra in — a solution to a problem in the domain 
of assembly planning. This application uses several types of operations on arrangements of 
geodesic arcs embedded on the sphere as basic blocks. In this context the chapter introduces 
exact implementations of additional applications that exploit geodesic arcs embedded on the 
sphere, such as polyhedra central-projection and Boolean set-operations applied to point sets 
embedded on the sphere bounded by geodesic arcs. Great parts of this chapter are extracts 
from a paper recently published in the 8^^ International Workshop on Algorithmic Founda- 
tions of Robotics (WAFR) |FH08] . Specific background of assembly planning is provided at 
the beginning of the chapter. 

We refer the reader to some ongoing research and future prospects and conclude in 
Chapter [6l 

The software access-information along with some further design details are provided in 
the Appendix. 
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A common mistake that peo- 
ple make when trying to 
design something completely 
foolproof is to underestimate 
the ingenuity of complete 
fools. 

Douglas Adams 




Arrangements on Surfaces 



Given a finite collection C of geometric objects (such as lines, planes, or spheres) the arrange- 
ment A{C) is the subdivision of the space where these objects reside into cells as induced 
by the objects in C. In this thesis we deal only with arrangements embedded on certain 
two-dimensional orientable parametric surfaces in M^, i.e., planes, cylinders, spheres, tori, 
and surfaces homeomorphic to them. In this case the objects in C embedded on the surface 
5* are curves that divide 5* into a finite number of cells of dimension (vertices), 1 (edges) 
and 2 (faces). Figure [2?1] shows various types of arrangements embedded on two-dimensional 
parametric surfaces. 

The Arrangement_on_surf ace_2^ package of Cgal |WFZH0i7a] is a generic implemen- 

^As a convention, Cgal prescribes the suffix _2 for all data structures of planar objects and the suffix _3 
for all data structures of 3D objects. In the case of arrangements on surfaces the suffix _2 indicates the 
dimension of the parameter space of the embedding surface. 




(a) (b) (c) 

Figure 2.1: (a) An arrangement of circles in the plane, (b) an arrangement of lines in the plane, and 
(c) an arrangement of geodesic arcs on the sphere. 
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tation of a complete software package that constructs and maintains arrangements em- 
bedded on two-dimensional parametric surfaces jFHK+07[ IFWH041 IWFZH07bl lBFH+07[ 
lBFH"'"09b] (see Section [1.3.31 for the history of the package). As arrangements are ubiq- 
uitous in computational geometry and have many theoretical and practical applications( 
see, e.g., |dBvKOS00l lASOOl IHal04] ) . many potential users in academia and in industry can 
benefit from the Arrangement_on_surf ace_2 package. 

As mentioned in Section ri.3.2I CGAL in general and the Arrangement_on_surf ace_2 pack- 
age of Cgal in particular follow the Exact Geometric Computation paradigm. The devel- 
oped code covers all cases to successfully handle degenerate input, while ideal computer 
arithmetic is emulated. While the Cgal package supports arrangements induced by general 
algebraic curves of arbitrary degree, the constructions and uses of arrangements described 
in this thesis require only rational arithmetic. This is a key property that enables efficient 
implementations of all the algorithms presented in the thesis. 

Cgal arrangements and their related components are the results of ongoing collaborative 
research we were, and still are, deeply involved with; see Section 11.3.31 for the evolution of 
the arrangement package. We, in particular, had a significant impact on specific topics 
within this research follows: We came up with great parts of the geometry-traits 

concept hierarchy, and the corresponding sets of minimal requirements (see Section 1^751) : we 
contributed the geometry-traits module for polylines (see Section I2.5.4p and the geometry 
and topology traits modules for geodesic arcs embedded on the sphere (see Section I^TBI) : we 
led the design of the Boolean set-operation architecture (See Section [2.7. ip . and we provided 
solutions to numerous issues encountered during the development of all these components. 
This chapter provides a comprehensive overview of these components, with a slight emphasis 
on asspects related to our work. It explains how arrangements induced by any type of curves 
can be constructed, maintained, and used by other applications. 



2.1 Related Work 

Both closely related Mapc [KCMKOOj and ESOLID |CKF+04] |H] libraries consist of 
an arrangement-construction module for algebraic curves. However, these implementations 
make some general-position assumptions on the input curves. The Leda library [MNOO] 
includes geometric facilities that allow the robust construction and maintenance of planar 
maps of line segments that may contain degeneracies. However, the resulting planar maps 
are represented as simple graphs that cannot fully describe the topological structure of the 
arrangement. For example, it is impossible to encode the containment relation between 
disconnected components of the graph (i.e., to keep track of the holes contained in a face; 
see Section [2. 3. 2p . LEDA-based implementation of arrangements of conic curves and of cubic 
curves were developed under the ExACUS project [BEH"'"05] [S]. 

Cgal's arrangement package was the first complete software implementation, designed 
for constructing arrangements of arbitrary planar curves and supporting operations and 
queries on such arrangements. The package was employed by many users to develop a 
variety of applications in various domains. For example, it was used to solve geometric op- 
timization problems |Rog03[ ICULYKT03j , to construct Minkowski sums efficiently jAFH02[ 



2.2. Parametric Surfaces 



19 



IFH071 IFSHOSb] , to design snap-sounding algorithms |HP02j , to construct envelopes of sur- 
faces |Mey06| . It was used in motion planning |HH03[rWvH07j . assembly partitioning [FH08], 
cartography [DHHOlj , and several other applications [H (TB] . The package was also used to 
compute arrangements of quadrics |BHK"'"05] by considering the planar arrangements of their 
projected intersection curves. A better approach to compute such arrangements |BFH+n7[ 
lBFH+n9aj was introduced once the package started to support arrangement on paramet- 
ric surfaces, which also enables the computation of arrangements embedded on Dupin cy- 
clides |BK08j . The torus, for example, is a Dupin cyclide. 

Sweeping the plane with a line is one of the most fundamental algorithmic mechanisms 
in computational geometry. The Arrangement_on_surf ace_2 package includes a generic 
implementation of an elaborate version of this mechanism exploited by several higher-level 
operations supported by the package. 

The famous sweep-line algorithm of Bentley and Ottmann [B079j was originally formu- 
lated for sets of non- vertical line segments, with the "general position" assumption that no 
three segments intersect at a common point and no two segments overlap. Many generaliza- 
tions have been introduced ever since, such as the ability to handle more general curves |SH89] 
and to deal with degeneracies (see |dBvKOS00l Section 2.1] and jMNOOl Section 10.7] for a 
discussion about degeneracies induced by line segments). 

Effective algorithms for manipulating arrangements of curves have been a topic of con- 
siderable interest in recent years with an emphasis on exactness and efficiency of implemen- 
tation |FHK+n7j . Mehlhorn and Seel |MSn 3; propose a general framework for extending the 
sweep-line algorithm to handle unbounded curves; however, their implementation can only 
handle lines in the plane. Arrangements on spheres are covered by Andrade and Stolfi [ASOlj . 
Halperin and Shelton |HS98j , and recently Cazals and Loriot |CL06j . Cazals and Loriot have 
developed a software package that can sweep over a sphere and compute exact arrangements 
of circles on it. 

The Leda external package (LEP) SphereGeometry [ISJ handles geodesic arcs on the 
sphere using an implicit representation, which enables the use of exact rational arithmetic 
to handle objects of this type. The package contains implementations of basic algorithms 
related to geodesic arcs on a sphere, such as computing the spherical convex hull, the union 
of two spherical polygons, and the width of a three-dimensional set of points. It does not, 
however, support arrangements. 

2.2 Parametric Surfaces 

A parametric surface S is defined by a continuous function : P — M^, where the domain 
F = U xV is a rectangular two-dimensional parameter space with bottom, top, left, and right 
boundaries, and the range S = /s'(P). U and V are open, half-open, or closed intervals with 
endpoints in MU {— oo, +oo}.We use Mmin, ^max, ^^min, and fmax to denote the endpoints of U 
and V, respectively. For example, the standard parameterization of the plane is fs{u,v) = 
{u,v,0), where U = V = (— cx), -t-cx)), and the unit sphere is commonly parameterized as 
fs{u, v) = (cosM cosf , sinw cosf , sinw), where P = [— vr, vr] x [— |, |]. 

A contraction point p G S* is a singular point, which is the mapping of a whole boundary of 
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the domain P. For example, if the top boundary is contracted, we have \/u G U, fs{u, fmax) = 
p' for some fixed point p' G M^. An identification curve C C 5* is a continuous curve, 
which is the mapping of opposite closed boundaries of the domain P. If the left and right 
boundaries are identified, we have Wv G V, /^(umin, f ) = fs{uraax,v), (and similarly for the 
bottom and top boundaries). For example, consider the sphere as parameterized above. Its 
contraction points are (0, 0, ±1), as fs{u, — f ) = (0, 0, —1) and fs{u, f ) = (0, 0, 1) for all u. 
Its identification curve is {/5(vr,t>) | — | < f < f}, as vr,f) = fs{+'^,v) for all v. 

A parameterizable curve 7 is a continuous function j : I —>■ F where / is an open, half- 
open, or closed interval with endpoints and 1, and 7 is injective, except for at a finite 
number of points. If ^ /, limt^o+ 7(^) exists (in the closure of P) and lies in an open 
side of the boundary. Similarly, if 1 ^ /, limj^i_ 7(t) exists and lies in an open side of the 
boundary. A curve C in S* is the image of a curve 7 in the domain. 

A curve is closed in the domain if 7(0) = 7(1); in particular, G / and 1 G /. A curve 
is closed in the surface S (or simply closed) if /s'(7(0)) = /5(7(1)). A curve 7 has two ends, 
the 0-end (7,0) and the 1-end (7, 1). If c? G /, the d-end has a geometric interpretation. It 
is a point in P. li d ^ I, the d-end has no geometric interpretation. You may think of it as a 
point on an open side of the domain or an initial or terminal segment of 7. If d ^ I, we say 
that the d-end of the curve is open. Consider for example the equator curve on the sphere 
as parameterized above. This curve is given by •yit) = (7r(2t — 1), 0), for t G [0, 1]. The 0-end 
of 7 is the point (— vr, 0) in P and a point on the equator of the sphere. It is closed on the 
sphere, but non-closed in P. 

A u-monotone curve is the image of a curve 7, such that if ti < t2, then M(7(ti)) < 
"^^(7(^2)) for ^1 < ^2- A vertical curve is the image of a curve 7, such that u{''j{t)) = c for all 
t G / and some c (zU and f (7(^1)) < f (7(^2)) for ti < ^2- For instance, every Meridian curve 
of a sphere parameterized as above is vertical. A weakly u-monotone curve is either vertical 
or M- monotone.^ 

The Arrangement_on_surf ace_2 package handles inducing curves that are decompos- 
able into parameterizable weakly w-monotone curves as defined above. Any two weakly 
M-monotone curves must intersect only a finite number of times or overlap only in a finite 
number of sections, if at all. The curves must be embedded on parameterizable surfaces 
as defined above. A curve can be unbounded or reach the boundaries of the embedding 
surface. A boundary may define a contraction point or an identification curve. ^ We allow 
non-injectivity on the boundary, denoted 9P, and require bijectivity only in P \ (9P (the 
interior of P). More precisely, we require that fs{ui,Vi) = fs{u2,V2) and {ui,Vi) 7^ {u2,V2) 
imply {ui,Vi) G dF and (^2,^2) G dF. Informally, we require that all geometric operations 
defined in Section 12.51 be applicable on our curves. 

Code reuse is maximized by generalizing the prevalent algorithms and their implementa- 
tions originally designed to operate on arrangements embedded in the plane. The generalized 
code handles features embedded in the modified surface S : fg = fs{u,v) | (M,f) G P \ dF 
defined over the interior of the parameter space, where identification curves, contraction 
points, and points at infinity are removed. Specific code that handles unbounded features 



^u-monotone curves refer to weakly u-monotone curves hereafter. 

■^We do not support surfaces that contain a contracted identification curve. 
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or features that reach the boundaries is added to yield a complete implementation. 

2.3 The Arrangement Package Architecture 

The main class of the package, namely Arrangement_on_surf ace_2, constructs and maintains 
the embedding of a set of continuous weakly u-monotone curves that are pairwise disjoint 
in their interiors on a two-dimensional parametric surface in R^. It provides the necessary 
capabilities for maintaining the embedded graph, while associating geometric data with the 
vertices, edges, and faces of the graph. The embedded graph is represented using a doubly- 
connected edge list (Dcel) [dBvKOS00[ Section 2.2], which maintains the incidence relations 
on its features |WFZH07b] . Each edge of the subdivision is represented by two halfedges 
with opposite orientation, and each halfedge is associated with the face to its left. It is 
based on an implementation of a halfedge data-structure (Hds) |Ket07b] also used by the 
polyhedral-surfaces package |Ket99l IKet07aj . 

An important guideline in the design is to decouple the arrangement representation from 
the various algorithms that operate on it. Thus, the Arrangement_on_surf ace_2 class pro- 
vides only a restricted set of methods for locally modifying the arrangement; see Section [2.3.2[ 
Non-trivial algorithms that involve geometric operations are implemented as free (global) 
functions that use the interface of the arrangement class; see Section 12. 4[ For example, the 
package offers free functions for incremental or aggregated insertion of curves that may not 
necessarily be w-monotone, and the insertion location of which are unknown a priori. Each 
input curve is subdivides into several w-monotone subcurves before inserted using one of the 
member methods listed in Section [2.3.21 

2.3.1 The Data Structure 

The Arrangement_on_surf ace_2<GeometryTraits ,TopologyTraits> class-template must 
be instantiated with two types as follows: 

• A geometry-traits class, which defines the abstract interface between the arrangement 
data-structure and the geometric primitives it uses. It is tailored to handle a specific 
family of curves, and it encapsulates implementation details, such as the number type 
used, the coordinate representation (i.e., Cartesian or homogeneous), the algebraic 
computation methods, and auxiliary data stored with the geometric objects, if present; 
see Section 12.51 for more details. 

• A topology-traits class, which adapts the underlying Dcel to the embedding modi- 
fied surface S. It determines whether the embedded surface is bounded, or otherwise 
whether a boundary defines a contraction point or an identification curve. If the induc- 
ing curves are confined to the modified parameter space, the tasks of the topology-traits 
class are minimal. However, in other cases it maintains the features that escape the 
modified parameter space P. 

The underlying DcEL in turn associates a point with each vertex and a w-monotone curve 
with each halfedge pair, where the geometric types of the point and the w-monotone curve 
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are defined by the geometry-traits class. Users may extend the default DcEL data-structure, 
in order to attach additional data to the DcEL records, as explained in Section I2.3.3I 

The Arrangement_on_surf ace_with_history_2 class-template represents an arrangement 
of general curves embedded on a two-dimensional parametric surface, and maintains the con- 
struction history of the arrangement. Input curves that induce the arrangement are split into 
ti-monotone subcurves that are pairwise disjoint in their interior, and these subcurves are 
the embeddings of the arrangement halfedges. While using the Arrangement_on_surf ace_2 
class we lose track of the connection between input curves and their final embeddings, in 
the Arrangement_on_surf ace_with_history_2 data-structure each embedded w-monotone 
curve is extended to store a pointer to the input curve associated with it, or a container of 
curve pointers in case the embedded w-monotone curve is an overlapping section of several 
input curves. 

The Arrangement_on_surf ace_with_history_2 class is a simple decorator'^ for 
Arrangement_on_surf ace_2. It inherits from an Arrangement_on_surf ace_2 class-template 
instantiated with a geometry-traits class that extends the w-monotone curve type. It also 
stores the set of input curves, and maintains a data structure that enables efficient traversal 
of all halfedges induced by given input curves. The cross-pointers between input curves 
and arrangement halfedges are maintained using an observer (see Section I2.4.4p that keeps 
track of each change that involves an arrangement halfedge and updates its underlying curve 
accordingly. 

Users can traverse the original curves of each arrangement halfedge, or iterate over all 
halfedges induced by a given input curve. Tracing back the curve (or curves) that induced 
an arrangement edge is essential in a variety of applications that use arrangements, such as 
robot motion planning; see, e.g., [HH03j . It is possible, for example, to efficiently remove a 
curve from the arrangement by deleting all edges it induces. 

Arrangements embedded in the plane are very common and, at least as far as the ar- 
rangement package of Cgal is concerned, have a longer history than their generalization 
for two-dimensional surfaces in M^. The Arrangement_2 class-template represents a planar 
subdivision. It maintains the embedding of continuous weakly -u-monotone curves in the 
xy plane, parameterized the natural way. That is, the two parameters u and v are directly 
mapped to x and y, respectively. Thus, w-monotonicity implies x-monotonicity and vice 
versa. The Arrangement_2 class-template is parameterized with a geometry-traits class and 
with a DCEL data-structure. It inherits from an Arrangement_on_surf ace_2 class-template 
instantiated with the geometry traits template parameter and with a specific topology-traits 
class suitable for the plane. The dedicated topology traits is instantiated with the DcEL 
template parameter. Similarly the Arrangement_with_history_2 class-template represents 
a planar subdivision, and maintains the construction history of the arrangement. 

The package offers various operations on arrangements stored in these representations, 
such as point location, insertion of curves, removal of curves, and overlay computation. 



"^A decorator attaches additional responsibilities to an object dynamically |GHJV95] . 
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The interface of Arrangement_on_surf ace_2 consists of various methods that enable the 
traversal of arrangement features. The class supplies iterators over its vertices, halfedges, or 
faces. The classes Vertex, Half edge, and Face, nested in the Arrangement_on_surf ace_2 
class, supply in turn methods for local traversals. For example, it is possible to visit all 
halfedges incident to a specific vertex. Halfedges stored in doubly-connected lists form 
chains. The chains define the inner and outer connected components of the boundary (CCB) 
of each face. A bounded face in the Arrangement_2 data structure has a single outer CCB 
representing the outer boundary of the face, and may have several inner CCBs representing 
its holes. However, a face in the general Arrangement_on_surf ace_2 data structure may 

Naturally, it is possible to traverse all 



have several inner and outer CCBs; see Section 12.6 
the halfedges along the inner and outer boundaries of a given face 

Alongside with the traversal methods, the arrangement 
class also supports several methods that modify the ar- 
rangement. The functions insert_in_f ace_interior(C,f ) 
(Figure [2.21 (a)). insert_f rom_lef t_vertex(C,v) (Figure [22] 
(b)), insert_f rom_right_vertex(C,v)) (Figure [2.21 (c)). and 
insert_at_vertices(C,vl ,v2) (Figure [521(d)) create an edge 
that corresponds to a w-monotone curve C, whose interior is 
disjoint from existing edges and vertices. The choice of which 
one to use depends on whether the curve endpoints are asso- 
ciated with existing non-isolated arrangement vertices: (i) If 
both curve endpoints do not correspond to any existing ver- 
tex, insert_in_f ace_interior is used to generate a new 
hole inside an existing face, (ii) If exactly one endpoint cor- 
responds to an existing DcEL vertex, one of the functions 
insert_f rom_lef t_vertex() or insert_f roni_right_vertex() 
is called, depending on which endpoint is associated with an 
existing vertex. It forms an "antenna" emanating from an ex- 
isting connected component, (iii) Otherwise, both endpoints 
correspond to existing vertices, and insert_at_vertices() is 
called to connect these vertices using a pair of twin halfedges. These functions hardly involve 
any geometric operations, if at all.^ They accept topologically related parameters, and use 
them to operate directly on the DcEL records, thus saving algebraic operations, which are 
especially expensive when high-degree curves are involved. 

Other modification methods included in the arrangement class enable users to split an 
edge into two, to merge two edges incident to a common vertex, and to remove an edge from 
the arrangement. It is also possible to insert a point in the interior of a given face, creating 
an isolated vertex that corresponds to this point, or to remove an isolated vertex from the 
arrangement. 




Figure 2.2: The arrangement 
immediate insertion methods. 
The newly inserted curve is 
drawn in bright blue. Vertices 
created as a result of the inser- 
tion are drawn in bright red. 



^Unless we force checking preconditions. In this case the precondition evaluation involves geometric 
computation. 
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2.3.3 Cell Extension 

As mentioned above the Arrangement_on_surf ace_2 is parameterized by a topological traits, 
which in turn is parameterized by a DcEL class. Users may extend the default DcEL data- 
structure, in order to attach additional data to the DcEL records. The default DcEL model 
simply associates a point with each DcEL vertex and a m- monotone curve with each halfedge 
pair. Although it is possible to store auxiliary data with the curves or points by extending 
their respective types (see Section [2.5.3p . it is sometimes necessary to extend the vertex, 
halfedge, or face topological features of the DcEL. Many times it is desired to associate 
extra data with the arrangement faces only. For example, when an arrangement represents 
the subdivision of a country into regions associated with their population density. In this 
case there is no alternative other than to extend the DcEL face, as there is no geometry-traits 
class entity that corresponds to an arrangement face. A similar mechanism for extending 
topological features with auxiliary attributes can be found in other components of Cgal, 
such as the triangulation packages |PY07j and the polyhedral-surfaces package |Ket07aj . 

2.4 The Arrangement Facilities 
2.4.1 Sweep Line 

The Arrangement_on_surf ace_2 package offers a generic implementation of the sweep-line 
algorithm [dBvKOSOOl Section 2.1] in form of a class template called Sweep_line_2. It 
handles any set of arbitrary w-monotone curves, and serves as the foundation of a family of 
concrete operations, such as computing all intersection points induced by a set of curves, 
constructing an arrangement of curves, aggregately inserting a set of curves into an existing 
arrangement, and computing the overlay of two arrangements. A concrete algorithm is 
realized through a sweep-line visitor, a template parameter of Sweep_line_2, which follows 
the visitor design-pattern [GHJV95] . and models the concept SweepLineVisitor_2. In this 
case, a visitor defines an operation based on the sweep-line algorithm to be performed on an 
arrangement without the need to change the arrangement structure. Users may introduce 
their own sweep based algorithms by implementing an appropriate visitor class. ^ 

Another parameter of the Sweep_line_2 class-template is the geometry-traits class, which 
must be instantiated with a model of the ArrangementXMonotoneTraits-2 concept; see Sec- 
tion 12.51 for the precise definition of this concept. It defines the minimal set of geometric 
primitives, among the other, required to perform the sweep-line algorithm briefly described 
next. 

An imaginary vertical curve is swept over the surface from left to right, transforming the 
static two-dimensional problem into a dynamic one- dimensional one. At each time during 
the sweep a subset of the input w-monotone curves intersect this vertical line in a certain 
order. The subset of curves and their order along the sweep line may change as the line 
moves along the w-axis, implying a change in the topology of the arrangement, only at a 



^The Boost Graph Library, for example, uses visitors |SLL021 Section 12.3] to support user-defined 
extensions to its fundamental graph algorithms. 
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finite number of event points, namely intersection points of two curves and left endpoints or 
riglit endpoints of arcs of curves. The event points, namely endpoints and all the intersec- 
tion points that have already been discovered, are stored in a wf-lexicographic order in a 
dynamic event queue, named the U -structure. The ordered sequence of segments intersect- 
ing the imaginary vertical line is stored in a dynamic structure called the V -structure. Both 
structures are maintained as balanced binary trees that enable their efficient maintenance 
using an advanced implementation of red-black trees |Wei05] . 

During the sweep-line process the event objects in the [/-structure are sorted lexico- 
graphically, and the subcurve objects are stored in the l^-structure in the same order as the 
lexicographic order of their intersection with the imaginary sweep-line. The Sweep_line_2 
class performs only the operations required to maintain the [/-structure and the V^-structure, 
while the visitor class is responsible for producing the actual output of the algorithm. When- 
ever the sweep-line class handles an event point p, it sends a notification to its visitor. Using 
this information, the visitor can access the subcurves incident to p and the neighboring 
subcurves from above and below. 



2.4.2 Map Overlay 

The map overlay of two planar subdivisions Si and S2 is a planar subdivision iS, such that 
there is a face / in iS if and only if there are faces /i and /2 in Si and S2 respectively, such 
that / is a maximal connected subset of /i fl /2 [dBvKOSOOl Section 2.3]. The overlay of 
two two-dimensional subdivisions embedded on a surface is defined similarly. 

The overlay of two given arrangements, conveniently referred to as the "blue" and the 
"red" arrangements, is straightforwardly implemented using a sweep-line visitor. A consol- 
idated set of the "blue" and "red" curves is processed, while the imaginary vertical line is 
swept over the surface. The -u-monotone curve type is extended with a color attribute (whose 
value is either BLUE or RED); see Section [2.5.31 Using the extended type we filter out un- 
necessary computations. For example, we ignore monochromatic intersections, and compute 
only red-blue intersection points (or overlaps). This way the arrangement of a consolidated 
set of "blue" and "red" curves is computed efficiently. 

The overlay visitor needs to construct a DcEL that properly represents the overlay of 
two input arrangements, the Dcel's of which are potentially independently extended (see 
Section r2.3.3p . A face in the overlay arrangement corresponds to overlapping regions of the 
blue and red faces. An edge in the overlay arrangement is due to a blue edge, a red edge, 
or an overlap of two differently colored edges. An overlay vertex is due to a blue vertex, 
a red vertex, a coincidence of two differently colored vertices, or an intersection of a blue 
and a red curve. In each case, the data associated with the overlay DcEL feature should 
be computed from the red and blue DcEL features that induce it. To this end, the overlay 
visitor is parameterized by an overlay-traits module, which defines the merging operations 
between various DcEL features, achieving maximum genericity and fiexibility for the users. 
The instantiated overlay traits models the Overlay Traits concept. The concept requires the 
provision of ten functions that handle all possible cases as follows: 

1. A new vertex v is induced by coinciding vertices Vr and v^. 
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2. A new vertex v is induced by a vertex Vr that lies on an edge Cb- 

3. An analogous case of a vertex Vb that lies on an edge e^. 

4. A new vertex v is induced by a vertex Vr that is contained in a face /h. 

5. An analogous case of a vertex Vb contained in a face fr- 

6. A new vertex v is induced by the intersection of two edges Cr and Cb- 

7. A new edge e is induced by the overlap of two edges Cr and Cb- 

8. A new edge e is induced by the an edge that is contained in a face fb- 

9. An analogous case of an edge Cb contained in a face fr- 

10. A new face / is induced by the overlap of two faces fr and fb- 

We apply the overlay operations in four different ways in this thesis; see Sections I3.2.2[ I3.3.2[ 
I5.3.6[ and 15.3.71 for the different applications. Each application requires the provision of a 
different set of the ten functions above. 

2.4.3 Zone Construction 

The zone |Hal04j of a w-monotone curve C in an arrangement A is the set of cells of A{C) 
intersected by the curve C. 

The Arrangement_on_surf ace_2 package includes the Arrangement_zone_2 class-template, 
which computes the zone of an arrangement. Similar to the Sweep_line_2 template, the 
Arrangement_zone_2 template is parameterized with a zone visitor, a model of the concept 
Zone Visitor_2, and it serves as the foundation of a family of concrete operations, such as in- 
serting a single curve into an arrangement and determining whether a query curve intersects 
with the curves of an arrangement. 

The zone of a curve C is computed by locating the left endpoint of C in the arrangement, 
and then "walking" along the curve towards the right endpoint, keeping track of the vertices, 
edges, and faces crossed on the way (see, for example, [dBvKOSOOl Section 8.3] for the 
computation of the zone of a line in an arrangement of lines). 

It is sometimes necessary to compute the zone of a curve in an arrangement without 
actually inserting the curve. In other situations, the entire zone is not required, as in the case 
of a process that only checks whether a query curve passes through an existing arrangement 
vertex; if the answer is positive, the process can terminate as soon as the vertex is located. 
While the sweep-line algorithm operates on a set of input w-monotone curves and its visitors 
can just use the notifications they receive to construct their output structures, the zone- 
computation algorithm operates on an arrangement object and its visitors may modify the 
same arrangement object as the computation progresses. This makes the interaction of the 
main class with its visitors slightly more intricate. 

2.4.4 Observers 

Some arrangement-based algorithms and applications should be bound to a specific arrange- 
ment instance and receive notifications on various topological changes this arrangement 
undergoes. This is not just a convenience, but crucial to the usability of the package, as it 
might be the only way for providing an algorithm with a certain input, such as data that 
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should be bound to the topological features of the arrangement, and is available only during 
construction; see Section 12.3.31 for an example. 

The Arrangement_on_surf ace_2 package supports a notification mechanism, which fol- 
lows the ofeeri'er design-pattern [ GHJV95j . In this case of one-to-many dependency a set of 
observes depend on a single arrangement, so that when the arrangement changes state, all 
its dependents are notified and updated automatically. Using this mechanism it is possible 
to attach any number of observer instances to a specific arrangement, such that all attached 
observers get notified on local and global changes the arrangement undergoes. 

The Arr_observer<Arrangement> class-template, parameterized by an arrangement type, 
stores a pointer to an arrangement object, and is capable of receiving notifications just be- 
fore a structural change occurs in the arrangement and immediately after such a change 
takes place. Hence, each notification comprises of a pair of "before" and "after" functions 
(e.g., bef ore_split_f aceO and af ter_split_f aceO). The Arr_observer class-template 
serves as a base class for other observer classes and defines a set of virtual notification func- 
tions, giving them all a default empty implementation. The interface of the base class is 
designed to capture all possible changes that arrangements can undergo, with a minimal set 
of topological events. 

The set of functions can be subdivided into three categories as follows: 

1. Notifiers of changes that affect the entire topological structure. Such changes occur 
when the arrangement is cleared or when it is assigned with the contents of another 
arrangement. 

2. Notifiers of a local change to the topological structure, such as the creation of a new 
vertex or an edge, the splitting of an edge or a face, the formation of a new hole inside 
a face, the removal of an edge, etc. 

3. Notifiers of a global change initiated by a free (global) function, and called by the 
free function (e.g., incremental or aggregate insert; see Section 12. 3p . This category 
consists of a single pair of notifiers, neither of them is called by methods of the 
Arrangement_on_surf ace_2 class-template itself. Issuing point-location queries (or 
any other queries for that matter) between the calls to the "before" and "after" func- 
tions of this pair is forbidden.'' 

See |WF05j for a detailed specification of the arrangement observer class sketched above. 

Each arrangement object stores a list of pointers to Arr_observer objects, and whenever 
one of the structural changes listed in the first two categories above is about to take place, 
the arrangement object invokes the appropriate function of each of its observers. It also does 
so immediately after the change has taken place. In addition, a free function may choose to 
trigger a similar notification, which falls under the third category above. 

^This constraint improves the efficiency of the maintenance of auxihary data structures for the relevant 
point-location strategies, which have to update their data structures according to the changes the arrange- 
ment undergoes (see Section 12.4.51 for more details). Since no point-location queries are issued between 
the invocation of bef ore_global_change () and af ter_global_change() , it is not necessary to perform an 
update each time a local topological change occurs, and it is possible to postpone the updates until after the 
global operation is completed. 
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In case the new observer is attached to a non-empty arrangement, its constructor may 
extract the relevant data from the non-empty arrangement using various traversal methods 
offered by the public interface of the Arrangement_on_surf ace_2 class, and update any 
internal data stored in the observer. This is necessary, for example, in case of the non- 
stateless point-location strategies, as shown in the next section. 

2.4.5 Point Location 

Point location is defined as follows: Given a point, find the arrangement cell that contains it. 
The Arrangement_on_surf ace_2 package provides the means to answer this query. Typically, 
the result of the point-location query is one of the arrangement faces, but in degenerate 
situations the query point can lie on an edge, or it may coincide with a vertex. Since 
the arrangement representation is decoupled from the algorithms that operate on it, the 
Arrangement_on_surf ace_2 class does not support point-location queries directly. Instead, 
the package provides a set of classes that are capable of answering such queries, all are models 
of the concept ArrangementPointLocation_2. Each model employs a different algorithm or 
strategy for answering queries. A model of this concept must define the locate () function, 
which accepts an input query point and returns an object representing the arrangement 
cell that contains this point (a polymorphic CGAL: : Object instance that can either be a 
Face_handle, a Half edge_handle, or a Vertex_handle). 

The following models of the concept ArrangementPointLocation_2 are included in the 
Arrangement_on_surf ace_2 package. 

• Arr_naive_point_location locates the query point naively, by exhaustively scanning 
all arrangement cells. It is the only strategy with unlimited support; see Section [2^61 

• Arr_walk_along_a_line_point_location simulates a reverse traversal along an imag- 
inary vertical ray emanating from the query point toward infinity. It starts from the 
unbounded face of the arrangement and moves downward toward the query point until 
it locates the arrangement cell containing it. 

• Arr_landmarks_point_location<Generator> uses an auxiliary generator class to cre- 
ate a set of "landmark" points, whose location in the arrangement is known. Given a 
query point, it uses a nearest-neighbor search structure (e.g., KD-tree) to find the near- 
est landmark, and then traverses the straight-line segment connecting this landmark 
to the query point.* See |HH08j for more details. 

• Arr_trapezoidal_ric_point_location implements Mulmuley's point-location algo- 
rithm |Mul90j . which is based on the vertical decomposition of the arrangement into 

^The "landmarks" strategy, requires that the arrangement is instantiated with a geometry-traits class 
that models the ArrangementLandmarksTraits-2 concept, which adds two requirements to the basic Arrange- 
mentBasicTraits-2 concept: (i) Approximating the coordinates of a given point p using the double-precision 
arithmetic, and (ii) constructing a w-monotone curve that connects two given points p and q, where p rep- 
resents a landmark point and q is the query point. Most traits classes included in the arrangement package 
are models of this refined concept. 



2.5. Geometry- Traits Concepts 



29 



pseudo-trapezoids, and maintains a history directed acyclic graph (DAG) on top of the 
decomposition. 

The last two strategies have query times that are shorter than the query times of the first two. 
However, they require preprocessing and consume more space, as they maintain auxiliary 
data structures. The first two strategies do not require any extra data and operate directly 
on the DcEL that represents the arrangement. For a complete survey see |HH08] . 

Each of the "landmarks" point-location class and the trapezoidal point-location class uses 
an observer to receive notifications whenever the arrangement is modified. For example, the 
default generator employed by the "landmarks" strategy uses the arrangement vertices as 
landmarks, so whenever a new vertex is created (by the insertion of a new edge, by the 
splitting of an existing edge, or by the insertion of an isolated point), it should be inserted 
into the nearest-neighbor search structure maintained by the respective landmark class. The 
usage of the notification mechanism makes it possible to associate several point-location 
objects with the same arrangement simultaneously. 

The "landmarks" and the trapezoidal point-location strategies are both characterized by 
very efficient query time at the cost of time-consuming preprocessing. Naturally, these strate- 
gies exhibit better overall performance when the number of arrangement updates is relatively 
small compared to the number of issued queries. For a report on extensive experiments with 
the various point-location strategies see [HH08j . 

2.5 Geometry- Traits Concepts 

The implementations of the 
various algorithms that con- 
struct and manipulate ar- 
rangements are generic, as 
they are independent of the 
type of curves they handle. 
All steps of the algorithms 
are enabled by the minimal 
set of geometric primitives 
gathered in the geometry- 
traits class, a model of a 
geometry-traits concept. The 
geometry-traits concept is 
factored into a hierarchy of 
refined concepts. The refine- 
ment hierarchy is defined according to the identified minimal requirements imposed by dif- 
ferent algorithms that operate on arrangements, thus alleviating the production of traits 
classes that handle complicated curves, and increasing the usability of the algorithms. The 
requirements listed by the geometry-traits concepts include only the utterly essential types 
and operations, and fully specify all the preconditions that the input must satisfy, as these 
may simplify the implementation of models of this concept even further. 



ArrangementBasic Traits-2 




ArrangementLandmarksBasicTraits-2 

Y 

ArrangementX Monotone Traits-2 




ArrangementLandmarksXMonotoneTraits-2 

Y 

ArrangementTraits-2 




ArrangementLandmarksTraits-2 



Figure 2.3: Refinement hierarchy of geometry-traits concepts. 
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The following sections are dedicated to a detailed description of the hierarchy. We list 
the minimal requirements of each layer in the hierarchy, and provide formal definitions for 
the required operations. The letters x and y are used in the code to refer to the two surface 
parameters, as arrangements embedded in the xy-plane are more common and familiar. 
The names of required nested types (e.g., X_monotone_curve_2) and valid expressions (e.g., 
compare_x) are faithful to the original source code. However, we use the letters u and v 
in the formal definitions below to refer to the two surface parameters, as these definitions 
apply to the general case of arrangements embedded on surfaces. Let cmp^() and cmp^() 
denote two predicates that accept two points and compare them by their w-coordinates and 
by their ^-coordinates respectively. We use the following notation. For a point p, {up,Vp), 
denotes a pre-image, and for a curve C, 7 denotes a pre-image, that is, p = fs{up,Vp) and 
C(t) = fsilit)) for all t e I. 

The basic concept ArrangementBasicTraits-2 requires the definition of the types Point_2 
and X_monotone_curve_2. The latter represents a w-monotone curve, and the former is the 
type of the endpoints of the curves, representing a point on the surface. This concept 
lists the minimal set of predicates on objects of these two types sufficient to enable the 
operations provided by the Arrangement_on_surf ace_2 class-template itself, namely the 
insertion of bounded w-monotone curves that are interior disjoint from any vertex and edge 
in the arrangement. All points and curves in the set below are required to have an inverse 
pre-image in P \ dF. In particular all curves are w-monotone. 



Compare_x_2: Compare two points by their u-coordinates. 

Compare_xy_2: Compare two points lexicographically by their u and then by their f-coordi- 
nates. 

Construct_min_2: Return the lexicographically smaller (left) endpoint of a given curve. 

Construct_max_2: Return the lexicographically larger (right) endpoint of a given curve. 

Is_vertical_2: Determine whether a weakly w-monotone curve is vertical. 

Compare_y_at_x: Given a point p and a curve C, such that the Up lies in the M-range of C, 
determine whether p is above, below, or lies on C. More precisely, if C is vertical, 
determine whether p lies on C, or above or below C. Otherwise, since ^(7(0)) <Up< 
14(7(1)) must hold and C is u-monotone, there must be a unique < t' < 1, that 
satisfies M(7(t')) = Up. Return cmp„(p, 7(t')). 

Compare_y_at_x_right: Given two curves Ci and C2 that share a common left endpoint p, 
determine the relative position of the two curves immediately to the right of p. More 
precisely, return cmp^ (71 (ei), 72(62)), where €1,62 > are infinitesimally small. 

Compare_y_at_x_lef t: Given two curves Ci and C2 that share a common right endpoint p, 
determine the relative position of the two curves immediately to the left of p. More 
precisely, return cmp^(7i(l — ei), 72(1 — £2)), where ei, €2 > are infinitesimally small. 
This is an optional requirement with ramifications in case it is not fulfilled; see Sec- 
tion ESU 



The set of predicates hsted above is also sufficient for answering point-location queries by the 
various point-location strategies, with a small exception of the "landmarks" strategy, which 
requires a traits class that models the refined concept ArrangementLandmarksTraits-2. This 
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is described in Section 12.4.51 

Constructing arrangements induced by w-monotone curves that may intersect in their 
interior, requires an arrangement instantiated with a traits class that models the concept 
ArrangementXMonotoneTraits-2 . This concept refines the basic arrangement-traits concept 
described above, as it requires an additional method for computing intersections between 
u-monotone curves, among the other. An intersection point between two curves is also rep- 
resented by the Point_2 type. The refined traits concept also requires methods for splitting 
curves at these intersection points to obtain pairs of interior disjoint subcurves and merging 
pairs of subcurves. In summary, a model of the refined concept must provide the additional 
operations bellow. All points and curves in the set below are required to have an inverse 
pre-image in P \ dW. In particular all curves are -u-monotone. 



Intersection_2: Compute the intersections between two given curves Ci and C2. 
Split_2: Split a given curve C at a given point p, which lies in the interior of C, into two 

interior disjoint subcurves. 
Merge_2: Merge two mergeable curves Ci and C2 into a single curve C. 
Is_mergeable_2: Determine whether two curves Ci and C2 that share a common endpoint 

can be merged into a single continuous curve representable by the traits class. 



The further refined concept ArrangementTraits-2 enables the construction of arrange- 
ments induced by general curves. A model of the refined concept must define a third type 
that represents a general (not necessarily w-monotone) curve, named Curve_2. It also has 
to supply a method that subdivides a given curve into simple u-monotone subcurves, and 
possibly isolated points.^ We refer to the entire hierarchy of refinements defined above as 
a single "abstract" concept called NoBoundary Traits, as it represents concepts the models 
of which handle curves that must have inverse pre-images in P \ dW. We use this abstract 
concept to simplify the description of the hierarchy defined below. 

The package introduces addi- 
tional concepts, models of which are 
able to handle unbounded curves or 
curves that reach the boundaries, the 
endpoints of which coincide with con- 
traction points or lie on identifica- 
tion curves; see Figure 12.41 The 
"abstract" Has Boundary Traits sub- 
hierarchy lists additional predicates 
required to handle both curves that 
reach or approach the boundaries of 
the parameter space. It has no mod- 
els. The refined BoundedBoundaryTraits and UnhoundedBoundary Traits sub-hierarchies list 
additional predicates required to handle bounded and unbounded curves, respectively. The 



^For example, the curve {x^ + y^)(a;^ + — 1) = is comprised of two u-monotone circular arcs, which 
together form the unit circle, and a singular isolated point at the origin. 



NoBoundary Traits 



HasBoundary Traits 



UnhoundedBoundary Traits BoundedBoundaryTraits 



AllBoundary Traits 



Figure 2.4: Abstract refinement hierarchy of geometry- 
traits concepts for arrangement on surfaces. 
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geometry-traits class that handles arcs of great circles models the BoundedBoundary Traits 
concept, as the parameter space is bounded in all four directions. Finally, the AllBoundary- 
Traits sub-hierarchy refines all the above. A model of this concept can handle unbounded 
curves in some directions and bounded curves in others. 

In the rest of this section all curves are required to be u-monotone. The HasBoundary- 
Traits concept requires the following additional operations: 



Parajneter_space_irL_x_2: Given a curve C and an index d G {0, 1} that identifies one of its 
ends, determine the location of its pre-image in the domain P along the u dimension. 
More precisely, determine whether u{'y{d)) is equal to Wmm, Wmax, or falls in between. 
In case of an unbounded curve, determine whether limt^du{'j{t)) is equal to —oo or 

+ 00. 

Parajneter_space_in_y_2: Given a curve C and an index d G {0, 1} that identifies one of its 
ends, determine the location of its pre-image in the domain P along the v dimension. 
More precisely, determine whether v{'j{d)) is equal to Wmin, "ymax, or falls in between. 
In case of an unbounded curve, determine whether limt^dvilit)) is equal to — oo or 

+ 00. 

Compare_x_near_boundary_2: There are two predicates: 

1. Given a point p, the inverse of which is in P \ dF, a curve C, and an index 
d E {0, 1} that identifies an end of C, compare the u coordinates of p and a point 
along C near its given end. More precisely, return cmp„(p, 7(|(i — e|)), where e > 
is infinitesimally small. 

2. Given two curves Ci and C2 and two corresponding indices di,d2 G {0,1} that 
identify two ends of Ci and C2 respectively, compare the u coordinates of two 
points along Ci and C2 respectively near their given ends. More precisely, return 
cmp„(7i(|(ii — ei I), 72(1^2 — £2!)), where ei,e2 > are infinitesimally small. 

See Section [2.6. II for an example. 
Compare_y_near_boundary_2: Given two curves Ci and C2, and a single index d G {0, 1} that 
identifies two ends of Ci and C2, compare the v coordinates of two points along Ci and 
C2 respectively near the given ends. More precisely, return cmp^(7i(|(i — ei|),72(|(i — 
€2!)), where €1,62 > are infinitesimally small. See Section [2.6. II for an example. 



The UnboundedBoundaryTraits concept requires the following additional operations: 



Is_bounded_2: Given a curve C and an index d G {0, 1} that identifies an end of C, determine 
whether the curve end is bounded. 



The BoundedBoundary Traits concept requires the following additional operations: 



Is_on_x_identif ication_2: This predicate applies only to a parameterization that has a 
vertical identification curve. Given a point p (respectively a curve C), determine 
whether p (respectively C) lies on the vertical and identified sides of the boundary. 
More precisely, determine whether Up G {wmin, ^^max}- (Respectively, determine whether 
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U(7(t)) e {Mmin,Mmax}, Vt G [0,1]. 

Is_on_y_identif ication_2: This predicate applies only to a parameterization that has a 
horizontal identification curve. Given a point p (respectively a curve C), determine 
whether p (respectively C) lies on the horizontal and identified sides of the boundary. 
More precisely, determine whether Vp G {fmin, "Wmax} for all pre-images of p. (Respec- 
tively, determine whether v{'y(t)) G {fmin, "Wmax}, Vt G [0, 1].) 

Is_on_x_contraction_2: This predicate applies only to a parameterization that has a con- 
tracted vertical boundary, determine whether p coincides with a contraction point. 
More precisely, determine whether Up is equal to Wmin or Wmax- 

Is_on_y_contraction_2: This predicate applies only to a parameterization that has a con- 
tracted horizontal boundary, determine whether p coincides with a contraction point. 
More precisely, determine whether Vp is equal to Wmin or fmax- 

Compare_x_on_identif ication_2: This predicate applies only to a parameterization that 
has a horizontal identified sides of the boundary. Given two points pi and p2 that lie 
on the horizontal identification arc, compare their w-coordinates. 

Compare_y_on_identif ication_2: This predicate applies only to a parameterization that 
has a vertical identified sides of the boundary.. Given two points pi and p2 that lie on 
the vertical identification arc, compare their f-coordinates. 



All traits-class operations are implemented as function objects (functors) according to 
Cgal's guidelines. This allows extending the geometric types above, without the need to 
redefine the methods that operate on them; see |HHK^07] for details on the extensible kernel. 
For a detailed specification of the various concept requirements see |WF05j . 



2.5.1 The Geometry- Traits Adaptor 

The geometry-traits adaptor class-template implements geometric operations that are not 
provided by a model of the geometry-traits concept itself, using the operations supplied by 
a model of the geometry-traits concept as basic blocks. It decreases the effort required to 
develop geometry-traits models, and at the same time increases the usability of the geometry- 
traits models, adapting them for extended uses. A geometry-traits type is injected as a 
template parameter into the adaptor class, which inherits from it, centralizing all geometric 
operations. In cases where the efficiency of methods is crucial, a developer has a way to 
override these methods with optimized ones. 

For example, in order to determine whether a point p is in the w-range of a w-monotone 
curve C, the adaptor simply compares p to the endpoints of C. It checks whether p lies to 
the right of the left endpoint and to the left of the right endpoint. 

In some cases, the geometry-traits adaptor class uses a tag- dispatching mechanism to 
select the appropriate implementation of a geometry-traits class operation. Tag dispatching 
is a technique that uses function overloading to dispatch a function at compile time, based 
on properties of the types of the arguments the function accepts [3]. This mechanism enables 
users to implement their traits class with a reduced or alternative set of operations. The 
adaptor respects the tags hsted below every geometry-traits class must define. 



34 



Chapter 2. Arrangements on Surfaces 



Has_lef t_category: A Boolean tag that indicates whether the traits class provides the 
predicate compare_y_at_x_lef t, which compares two u-monotone curves to the left 
of a common right endpoint. This predicate is required only by some point-location 
strategies and by the zone-computation algorithm. While in some cases it is fairly easy 
for the traits-class implementer to provide it, in other cases it can be rather difficult, 
or even quite impossible. When this tag is false, the traits-class adaptor resorts to a 
somewhat less efficient algorithm that uses (other) existing traits-class predicates. 

Has_merge_category A Boolean tag that indicates whether a model of the ArrangementX- 
MonotoneTraits-2 supports the merge of w-monotone curves. If the tag is true, the 
traits class must provide the two operations merge_2 and is_mergeable_2. The merger 
operation is used to eliminate redundant features in the arrangement. For example, 
if we have a T-shaped structure formed by two line segments, and the vertical seg- 
ment forming the "leg" is removed, then it is possible to merge the two horizontal 
sub-segments. When the has-merge tag is false, the adaptor simply declares any pair 
of curves as non-mergeable. The only effect on the arrangement is that we cannot 
remove redundant vertices (of degree two) following the deletion of edges. 

Boundary_category: A quadruple tag that categorizes the traits class according to the hi- 
erarchy described in Figure 12. 4[ The adaptor provides empty implementations of the 
operations that are never invoked, yet required for smooth compilation. 

2.5.2 Geometry- Traits Models 



Table 2.1: Geometry-traits models 



Curve Family 


Degree 


Surface 


Boundness 


Arithmetic 


linear segment 


1 


plane 


bounded 


rational 


linear segments, rays, lines 


1 


plane 


unbounded 


rational 


piecewise linear curves 


oo 


plane 


bounded 


rational 


circular arcs, linear segments 


< 2 


plane 


bounded 


rational 


algebraic curves 


< 2 


plane 


unbounded 


algebraic 


quadric projections 


< 4 


plane 


unbounded 


algebraic 


algebraic curves 


< 3 


plane 


unbounded 


algebraic 


algebraic curves 


< n 


plane 


unbounded 


algebraic 


planar Bezier curves 


< n 


plane 


unbounded 


algebraic 


univariate polynomials 


< n 


plane 


unbounded 


algebraic 


geodesic arcs on sphere 


< 2 


sphere 


bounded 


rational 


quadric intersection arcs 


< 4 


quadric 


unbounded 


algebraic 


Dupin cyclide intersection arcs 


< n 


Dupin cyclides 


bounded 


algebraic 



The large number of geometry-traits models already implemented enables the construc- 
tion and maintenance of arrangements induced by many different types of curves. The 
package itself contains several models of the geometry-traits concept. A few other models 
have been developed by other groups of researchers. Models are distinguished not only by 
the different families of curve they handle, but also by their suitability for constructing and 
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maintaining arrangements with different characteristics. For example, there are two distinct 
models that handle line segments [WFZHOTbj . One caches information in the curve records, 
while the other retains the minimal amount of data. While operations on arrangements 
instantiated with the former model consume more space, they are more efficient for dense 
arrangements (namely, arrangements induced by curves with a large number of intersections). 
Another model handles not only (bounded) line-segments, but also rays and lines |BFH+07[ 
lBFH+09bj . There are traits models for non-linear curves, such as circular arcs [dCPTOTj . 
conic curves |Wei02t lBEH+02[ lEKP+04 ]. cubic curves |EKSW04j . and quartic curves that 
are the projection of the intersection of two quadric surfaces |BHK"'"05] . and there are traits 
classes for arcs of graphs of rational univariate polynomial functions |LPT08| IWFZHOTb] . 
There is even a traits class that handles algebraic curves of arbitrary degrees |EK08j . There 
is also a traits class that handles Bezier curves |HW07] . There is a traits class for geodesic 
arcs embedded on the sphere |FSH08b] lBFH"'"09aj . (see Section 12.61) . another one for in- 
tersections of quadrics embedded on a quadric |BFH+07[ lBFH+09aj . and another one for 
intersections of arbitrary algebraic surfaces with a Dupin cyclide embedded on the Dupin 
cyclide |BK081 lBFH"'"09a] . Finally, there is a model that handles continuous piecewise linear 
curves, referred to as polylines, (see Section [2.5.4p . 



2.5.3 Geometry- Traits Extension 

Traits-class decorators are used to extend the geometric entities defined by the traits class 
with additional, possibly non-geometric, data. An alternative way to achieve this is to extend 
the geometric types of the kernel, as the kernel is fully adaptable and extensible [HHK"'"07] . 
However, this indiscriminating extension may lead to an undue space-consumption, as every 
geometric object is extended, regardless of its use. It also requires nontrivial knowledge 
about the kernel structure and the techniques to extend it. 

There is a decorator that enables the extension of the (general) curve and the w-monotone 
curve types with distinct types of data, and there is a convenient one, where the data attached 
to the M-monotone curve type is a set of objects, the type of which is attached to the (general) 
curve type. This set usually contains a single data object, unless the w-monotone curve 
corresponds to an overlapping section of two curves or more. When a curve with a data field 
d is split into w-monotone subcurves, each subcurve is associated with a singleton set {d}. 
When two -u-monotone curves overlap, the decorator takes the union of their data sets, and 
associates it with the resulting overlapping subcurve. 



2.5.4 A Geometry-Traits Model that Handles Polylines 

Polylines are of particular interest, as they can be used to approximate more complex curves. 
At the same time handling them is easier than handling higher-degree algebraic curves, as 
rational arithmetic is sufficient to carry out exact computations on polylines. 

The geometry-traits model that handles polylines is a class-template called 
Arr_polyline_traits_2. It must be instantiated with a geometry-traits class that models 
the concept ArrangementLinear Traits. This concept refines the ArrangementTraits concept, 
as it adds a variant — it must handle line segments. This variant cannot be enforced by the 
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compiler, but rather be verified at run time. A polyline curve is represented as a vector of 
SegmentTraits : :X_monotone_curve_2 objects (namely segments). The polyline-traits class 
does not perform any geometric operations directly. Instead, it solely relies on the function- 
ality of the instantiated segment-traits class. For example, when we need to determine the 
position of a point with respect to a w-monotone polyline, we use binary search to locate 
the relevant segment that contains the point in its w-range, then we compute the position of 
the point with respect to this segment. Thus, operations on -u-monotone polylines of size m 
typically take O(logm) time. 

Users are free to choose the underlying segment-traits class based on the number of 
expected intersection points (see discussion above in Section [2.5.2p . Moreover, it is possible 
to instantiate the polyline-traits class-template with a traits class that handles segments 
with some additional data attached to them (see Section [2.5.31) . This makes it possible to 
associate different data objects with the different segments that compose a polyline. 



2.6 Arrangements of Geodesic Arcs on the Sphere 

In this section we concentrate on the particular category of arrangements embedded on 
surfaces, where the embedding space is the sphere, and the inducing objects are geodesic arcs. 
There is an analogy between this class of arrangements and the class of planar arrangements 
induced by linear curves (i.e., segments, rays, and lines), as properties of linear curves in the 
plane can be often, but not always, adapted to geodesic arcs on the sphere. 

An arrangement of geodesic arcs embedded on the sphere is defined as an instance of 
the Arrangement_on_surf ace_2 class-template instantiated with appropriate geometry- and 
topology-traits classes, namely Arr_geodesic_arc_on_sphere_traits_2 and 
Arr_spherical_topology_traits_2, respectively. The geometry-traits class is tailored to 
handle geodesic circs cls efficiently as possible, and defines the parameterization used: P = 
[—71 + a, IT -|- a] X [— f,f], fs{u,v) = (cosm cosw, sinw cosw, sint>), where a is a variable 
that must be set at compile time, and is by default 0. This parameterization induces two 
contraction points Ps = (0,0,-1) and Pn = (0,0,1), referred to as the south and north 
pole respectively, and an identification curve that coincides with the opposite Prime (Green- 
wich) Meridian. We developed the topology-traits class to support not only arrangements 
of geodesic arcs, but any type of curves embedded on the sphere parameterized as above, 
without compromising the performance of the operations gathered in the traits class. We 
hope that this topology-traits class will come in handy in the future for constructing and 
maintaining arrangements induced by types other than geodesic arcs, such as general cir- 
cular arcs, which appear in arrangements induced by intersections of spheres embedded on 
the sphere |CL06] . The topology-traits class initializes the DcEL to have a single face, the 
embedding of which, is the entire sphere. It is designed to retain the variant that this face 
always contains the north pole during modifications the arrangement may undergo. The 
topology-traits class is required, for example, to inform its users that the top and bottom 
boundaries of the parameter space are contracted and the left and right boundaries are iden- 
tified. It maintains a search structure of vertices that coincide with the contraction points 
or lie on the identification arc. 
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The figure to the right is a snapshot of a movie |FSH08b] that 
demonstrates, among the other, the sweep-line procedure carried out 
on the sphere. The red vertical arc that connects the poles is the 
imaginary sweep-line. The yellow vertex have been processed already. 
The dark blue arcs are curves that have been processed already and 
inserted into the arrangement. The light blue arcs are curves that are 
to be processed as the sweep line advances. 

2.6.1 The Geometry-Traits Model 

The geometry-traits class for geodesic arcs on the sphere is parameterized with a geometric 
kernel |HHK+07j that encapsulates the number type used to represent coordinates of geo- 
metric objects and to carry out algebraic operations on those objects. The implementation 
handles all degeneracies, and is exact, as long as the underlying number type is rational, 
even though the embedding surface is a sphere. We are able to use high-performance ker- 
nel models instantiated with exact rational number-types for the implementation of this 
geometry-traits class, as exact rational arithmetic suffices to carry out all necessary alge- 
braic operations. The ability to robustly construct arrangements of geodesic arcs on the 
sphere, and robustly apply operations on them using only (exact) rational arithmetic is a 
key property that enables an efficient implementation. 

A point in our arrangement is defined to be an unnormalized vector in M^, representing 
the place where the ray emanating from the origin in the relevant direction pierces the sphere. 
An arc of a great circle is represented by its two endpoints, and by the plane that contains 
the endpoint vectord and goes through the origin. The orientation of the plane and the 
source and target points determine which one of the two great arcs is considered. 

The point type is extended with an enumeration that indicates whether the vector (i) 
pierces the south pole, (ii) pierces the north pole, (iii) intersects the identification arc, or (iv) 
is in any other direction. An arc of a great circle is extended with three Boolean fiags that 
indicate whether any one of the x, y, z coordinates of the normal of the plane that defines the 
arc vanishes. These fiags are used to minimize the number of invocations of the geometry- 
traits operations, which has a drastic effect on the performance of arrangement operations 
at the account of a slight increase in space consumption. This representation enables an 
exact yet efficient implementation of all geometric operations required by the geometry-traits 
concept using exact rational arithmetic, as normalizing vectors (that represent directions and 
plane normals) is completely avoided. 

We describe in details four predicates, namely Compare_x_2, Compare_xy_2, 
Compare_x_near_boundary_2, and Compare_y_near_boundary_2; see Section [231 for the com- 
plete set of the concept requirements. The former compares two points pi and p2 by their 
M-coordinates. The concept admits the assumption that the input points do not coincide 
with the contraction points and do not lie on the identification arc. Recall that points are 
in fact unnormalized vectors that represent directions in R'^. We project pi and p2 onto 
the xy-plane to obtain two-dimensional unnormalized vectors pi and p2, respectively. We 
compute the intersection between the identification arc and the xy-plane to obtain a third 
two-dimensional unnormalized vector d. Finally, we test whether d is reached strictly before 
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P2 is reached, while rotating counterclockwise starting at pi. This 
geometric operation is supported by every geometric kernel of Cgal. 
In the figure on the right d is reached strictly before p2 is reached. 
Therefore, the w-coordinate of pi is larger than the w-coordinate of 

The predicate Compare_xy_2 compares two points pi and p2 lexico- 
graphically. It first applies Compare_x_2 to compare the u-coordinates 
of the two points. If the u-coordinates are equal, it applies a pred- ^ 
icate that compares the ^-coordinates of two points with identical 
M-coordinates, referred to as Compare_y_2. This predicate first com- 
pares the signs of the 2;-coordinates of the two unnormalized input ^ 
vectors. If they differ, it concludes that the point with the positive 
2;-coordinate has a f-coordinate that is larger than the f-coordinate 
of the point with the negative ^-coordinate. If the signs are identical, 
it compares the squares of their normalized z-coordinates, essentially 
avoiding the square-root operation. If the sign of the (unnormalized) z-coordinates of both 
points is positive (resp. negative), the point with the larger (resp. smaller) square of nor- 
malized z-coordinate has a larger f -coordinate. 

The predicates above accept points, the pre-images of which, lie 
in the interior of the parameter space. However, there is also a 
need to lexicographically compare the ends of arcs, the pre-images 
of which reach the boundary of the parameter space. The predicate 
Compare_x_near_boundary_2 accepts either (i) a point, the pre-image 
of which lies in the interior of the parameter space, and an arc end, 
or (ii) two arc ends. Such an arc end is provided by a vertical arc 
and an index that identifies one of the two ends of the arc, and must 
coincide with one of the contraction points. The first variant compares the w-coordinates 
of the input point and a point along the input arc near its given end, whereas the second 
variant compares the M-coordinates of two points along the input arcs near their respective 
given ends. Recall, that the w-coordinates of all points along a vertical arc are the same (C4 
and C5 in the figure above). Thus, we can compare the w-coordinate of an arbitrary point 
on a vertical arc that lie inside the parameter space. For example, for the second case, we 
compare the two vectors perpendicular to the normals to the planes that define the vertical 
arcs, respectively, e.g., the u-coordinate of a point on the arc C4 near its top end is smaller 
than the w-coordinate of a point on the arc C5 near its top end, and in particular it is smaller 
than the M-coordinate of the bottom end of C5. The Compare_y_near_boundary_2 predicate 
compares the f-coordinate of two arcs ends, the pre-images of which lie on the same (left or 
right) identified side of the boundary of the parameter space. We use the aforementioned 
Compare_y_2 predicate to compare the end points. If the points are equal, we compare the 
normals to the plane that define the arcs. In our example, the left end of Ci is smaller than 
the left end of C2, which is smaller than the left end of C3. 

All the required geometric types listed in the traits concept are maintained using only 
rational numbers. All required geometric operations are implemented using only rational 
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arithmetic.^" Degeneracies, such as overlapping arcs that occur during intersection compu- 
tation, are properly handled. The end result is a robust yet efficient implementation. 



2.7 Applications 

Arrangement on surfaces have many applications this thesis falls short to list. However, 
we do list a few samples we were involved (or remotely involved) with the implementation 
of which, i.e.. Regularized Boolean Set-Operations, Envelopes of Surfaces, and Voronoi di- 
agrams. Minkowski sum construction is covered in details in the following chapter. The 
Boolean set-operation results, minimization diagrams, maximization diagrams and Voronoi 
diagrams, (see Section [2.7.21 for definitions), and Minkowski-sums are all represented as ar- 
rangements, and as such can be passed as input to consecutive operations on arrangements 
supported by the Arrangement_on_surf ace_2 package and its derivatives. 



2.7.1 Regularized Boolean Set-Operations 



Together with R. Wein and B. Zuckerman we have developed a package that supports Boolean 
set-operations on point sets bounded by w-monotone curves embedded on two-dimensional 
parametric surfaces in |FWZH07j . In particular, it contains the implementation of reg- 
ularized Boolean set-operations, intersection predicates, and point containment predicates. 
A regularized Boolean set-operation op* can be obtained by ffist taking the interior of the 
resulting point-set of an ordinary Boolean set-operation (P op Q) and then by taking the 
closure |Hof04] . That is, P op* Q = closure (interior (P op Q)). Regularized Boolean set- 
operations appear in constructive solid geometry (CSG), because regular sets are closed under 
regularized Boolean set-operations, and because regularization eliminates lower dimensional 
features, namely isolated vertices and "antennas" (namely, dangling edges), thus simplifying 
and restricting the representation to physically meaningful solids. Ordinary Boolean set- 
operations, which distinguish between the interior and the boundary of a polygon, are not 
implemented within this package. However, we implemented a specialized ordinary union 
operation as part of an assembly partitioning application; see Chapter [51 

The operands and results of the regularized 
operations are general polygons that may have 
holes. The boundaries of a general polygon 
and of holes, if present, are general u-monotone 
curves. The Arrangement_on_surf ace_2 class is 
employed to represent a point set embedded on 
a two-dimensional parametric surface as an ar- 
rangement. A point set is typically constructed 
from a single general polygon or a collection of 
interior disjoint general polygons. The underly- 
ing arrangement must be instantiated with a geometry traits that models the concept Gener- 



ArrangementBasic Traits-2 



I 



ArrangementDirectionalXMonotoneTraits-2 



GeneralPolygonSetTraits_2 



Figure 2.5: Refinement hierarchy of geometry 
traits concepts for Boolean set-operations. 



""^^Points are represented as unnormalized vectors; The coordinates of such points are converted into ma- 
chine floating-point only for rendering purposes. 



40 



Chapter 2. Arrangements on Surfaces 



alPolygonSetTraits-2. This concept refines the concept ArrangementDirectionalXMonotone- 
Traits-2, which in turns refines the concept Arrangements asicTraits-2 (see Section [23]) . 

The ArrangementDirectionalXMonotoneTraits-2 concept treats its m- monotone curves 
as objects directed from one endpoint appointed to be the source to the other endpoint 
appointed to be the target. Thus, it requires few additional operations on w-monotone 
curves: 



Compare_endpoints_xy_2: Given a w-monotone curve C, compare the source and target 

points of C lexicographically. 
Construct_opposite_2: Given a u- monotone curve C, construct the opposite curve of C 

(namely, swap the source and target endpoints of C). 
Intersection_2: Compute the intersections between two given curves Ci and C2. 
Merge_2: Merge two mergeable curves Ci and C2 into a single curve C. 
Is_mergeable_2: Determine whether two curves Ci and C2 that share a common endpoint 

can be merged into a single continuous curve representable by the traits class. 



Most traits classes bundled in the Arrangement_on_surf ace_2 package and distributed with 
Cgal, are models of the concept ArrangementDirectionalXMonotoneTraits-2}^ 

The GeneralPolygonSetTraits-2 concept requires its models to define a type that rep- 
resents a general polygon and another one that represents general polygon with holes in 
addition to the Point_2 and X_monotone_curve_2 types that must be defined by all models 
of the generalized concept. It also requires the provision of several operations that operate 
on these two types listed below. 



Construct_polygon_2: Given a sequence C of w-monotone curves, construct a general poly- 
gon that has C as its outer boundary. 

Construct_curves_2: Given a general polygon P, obtain the sequence of w-monotone curves 
that comprise the boundary of P. 

Construct_general_polygon_with_holes_2: Given a general polygon P and a (possibly 
empty) set of holes 7i, construct a general polygon with holes that has P as its outer 
boundary and H, as its holes. 

Construct_outer_boundary: Given a general polygon-with-holes P, obtain the general poly- 
gon that is its outer boundary. 

Construct_holes: Given a general polygon-with-holes P, obtain the holes of P if any. 

Is_unbounded: Given a general polygon-with-holes P, determine whether it has an outer 
boundary. 



^^The Arr_polyline_traits_2 traits class is not a model of the ArrangementDirectionalXMonotone- 
Traits_2 concept, as the u-monotone curve it defines is always directed from left to right. Thus, an opposite 
curve cannot be constructed. 
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2.7.2 Envelopes 

Lower envelopes of functions on parametric surfaces are defined in a way similar to the 
standard definition of lower envelopes of bivariate functions in space |Hal04] . Let 5* be a 
two-dimensional parametric surface is M^. Given a set of bivariate functions F = {/i, . . . , /„}, 
where /j : S* — * M, their lower envelope \E'(M,f) is defined to be their point-wise minimum 
— mini<j<„ f). The minimization diagram J^[F) of the set F is the two- 
dimensional map obtained by the projection of the lower envelope onto S. Upper envelopes 
and maximization diagrams are defined analogously for the point-wise maximum of the 
functions. 




(a) (b) (c) 

Figure 2.6: Lower envelopes of various types of surfaces (a) The lower envelope of triangles, (b) The 
lower envelope of spheres, (c) The lower envelope of planes. (The minimization diagrams is drawn 
above the planes for clarity.) 



The Envelope_3 package of Cgal |Mey06 , IMWZ07] computes the lower (or the upper) 
envelope of a set of surfaces in R^. It is based on the Arrangement_on_surf ace_2 package, 
and like the base package, it handles degenerate input, and produces exact results. An ar- 
rangement data-structure is used to represent the resulting minimization diagram |Hal04j . 
The envelope computation is enabled by a traits class — a model of the concept Envelope- 
Traits-3, which refines the ArrangementTraits-2 concept. The Envelope_3 package currently 
contains three models of the EnvelopeTraits-3 concept that can be used to compute the en- 
velope of triangles, spheres, and planes, respectively; see Figure [2T6] for an illustration. Other 
models of the EnvelopeTraitsS concept have been developed, for example, a traits class that 
enables the computation of the envelope of a set of quadric surfaces |BM07j . 



2.7.3 Voronoi Diagrams 

Voronoi diagrams were thoroughly investigated, and were used to solve many geometric 
problems, since introduced by Shamos and Hoey to the field of computer science |SH75j 
(although their origin dates back centuries ago; see [UBSCOO] ). The concept of computing 
cells of points that are closer to a certain object than to any other object, among a finite 
number of objects, was extended to various kinds of geometric sites, ambient spaces, and 
distance functions, e.g., power diagrams of circles in the plane, multiplicatively weighted 
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Voronoi diagrams, additively weighted Voronoi diagrams |AK00j . This space decomposition 
is strongly connected to arrangements |ES86j . a property that yields a very general approach 
for computing Voronoi diagrams. 

Given a set of n points P = {pi, . . . ,Pn}, Pi £ S, we define R{P,Pi) = {x E S \ p{x,pi) < 
p{x,pj),j 7^ i}, where p{pi,Pj) is some given distance function. R{P,pi) is the region 
of all points that are closer to pi then to any other point in P. The Voronoi diagram of 
P over S is defined to be the regions R{P,pi), R{P,p2), ■ ■ ■ , R{P,Pn) and their boundaries. 
Edelsbrunner and Seidel j ES86j observed the connection between Voronoi diagrams in 
and lower envelopes in W^^^ of the corresponding distance functions to the sites. From the 
above definitions it is clear that ii fi : S ^ M. is set to be fi{x) = p{x,pi), for i = 1, . . . , n, 
then the minimization diagram of {/i, . . . , /„} over S is exactly the Voronoi diagram of P 
over S. 

K. E. Hoff et al. developed a technique for 
computing generalized 2D and 3D Voronoi dia- 
grams using interpolation-based polygon raster- 
ization hardware |HKL+99j . The hardware is 
used to draw the discrete and approximate lower 
envelope of the site-distance functions. Follow- 
ing similar principles, O. Setter et al. developed 
a new framework to compute different types of 
Voronoi diagrams embedded on certain paramet- 
ric surfaces in an exact manner [SSHOSj . The 
framework is based on the exact computation of 
the lower envelope of the site-distance functions 
over the surface |Mey06| . It provides a reduced 
and convenient interface between the construction of the diagrams and the construction of 
envelopes, which in turn are computed using the Envelope_3 package [MWZOTj . Obtaining 
a new type of Voronoi diagrams only amounts to the provision of a geometry-traits class 
that handles the type of bisector curves of the new diagram type. Essentially, every type of 
Voronoi diagram, the bisectors of which can be handled by a geometry traits class, can be 
implemented using this framework. In particular the geometry-traits class for geodesic arcs 
embedded on the sphere enables Voronoi diagrams of points on the sphere and their gener- 
alization, power diagrams, also known as Laguerre Voronoi diagrams, on the sphere, as the 
bisector curves between point sites on the sphere are great circles |NLC02l lOBSCOO] , and so 
are the bisectors between circle sites on the sphere under the Laguerre distance |Sug02] ; see 
Figure 12.71 Power diagrams on the sphere have several applications similar to the applica- 
tions of power diagrams in the plane. For example, determining whether a point is included 
in the union of circles on the sphere, and finding the boundary of the union of circles on the 
sphere tLLMSSj [5ug02l . 



(a) (b) 

Figure 2.7: Voronoi diagrams on the sphere. 

(a) The Voronoi diagram of 14 random points. 

(b) The power diagram of 10 random circles. 



^^In certain cases, the distance to a site may depend on various parameters associated with the site, e.g., 
in the cases of Mobius diagrams or anisotropic diagrams. 
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We implicitly construct envelopes of distance functions defined 
over the sphere to compute Voronoi diagrams. The image to the 
right illustrates the distance function from (0, 0) G [— vr, tt] x [— |, |] 
on the sphere in the parameter space. The great circle bisector of 
two point sites on the sphere is the intersection of the sphere and 
the bisector plane of the points in (imposed by the Euclidean 
metric). 

If a point on the sphere is given as a general vector in M^, it must be normalized. If 
the normalization results with a point with irrational coordinates, then it must be approx- 
imated to a point that lies exactly on the sphere |CDR92] . Once approximated, all further 
computations are carried out using exact rational arithmetic. 

Figure 12.81 (a) on the left shows an arrange- 
ment on the sphere induced by (i) the conti- 
nents and some of the islands on earth, and (ii) 
the institutions that hosted SoCG during this 
millennium, which appear as isolated vertices. 
The sphere is oriented such that College Park, 
MD, USA is at the center. The arrangement 
consists of 1054 vertices, 1081 edges, and 117 
faces. The data was taken from gnuplot [T3j and 
google maps [14J. Figure 12.81 (b) shows an ar- 
rangement that represents the Voronoi diagram 
of the nine cities, the institutions above are lo- 
Sedona, Pisa, New York, San Diego, Barcelona, 




(a) (b) 
Figure 2.8: Arrangements on the sphere. 



cated at, namely College Park, Gyeongju 
Medford, and Hong Kong. 

As mention above Voronoi diagrams, among the other, are repre- 
sented as arrangements and can be passed as input to consecutive oper- 
ations on arrangements supported by the Arrangement_on_surf ace_2 
package and its derivatives. The figure on the right shows the overlay 
of the two arrangements shown in Figure 12. 8[ 
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Efficiency is just intelligent 
laziness. 

Anonymous 




Minkowski Sum Construction 



We present two exact and robust implementations of efficient output-sensitive algorithms to 
compute the Minkowski sum of two polytopes in M^. We demonstrate the effectiveness of our 
Minkowski-sum computations through simple applications that exploit these operations to 
detect collision, and compute the Euclidean separation distance between, and the directional 
penetration depth of, two polytopes in M^. In Chapter Owe show a more involved application 
of these operations. 

Each method we have developed uses a different variant of Gaussian maps to main- 
tain dual representations of polytopes. Each method employs a different variant of two- 
dimensional arrangements to maintain the dual representations, and it makes use of many 
operations applied to arrangements in the corresponding representations. The first method 
uses the traditional (spherical) Gaussian map. The map is represented as an arrangement 
of geodesic arcs embedded on the unit sphere. The second method uses a data structure 
called Cubical Gaussian Map . It consists of six arrangements induced by linear segments 
embedded in the plane. The six arrangements correspond to the six faces of the unit cube 
— the parallel-axis cube circumscribing the unit sphere. 

A simple method to compute the Minkowski sum of two polytopes is to compute the con- 
vex hull of the pairwise sum of the vertices of the two polytopes. Although there are many 
implementations of various algorithms to compute Minkowski sums and answer proximity 
queries, we are unaware of the existence of complete implementations of methods to compute 
exact Minkowski sums other than (i) the naive method above, (ii) a method based on Nef 
polyhedra embedded on the sphere |HKM07] . and (iii) an implementation by Weibel [22] of 
Fukuda's algorithm |Fuk04] . Both our methods exhibit much better performance than the 
other methods in all demonstrated by the experiments reported in Table 13.51 Our 

methods well handle degenerate cases that require special treatment when alternative rep- 
resentations are used. For example, the case of two parallel facets facing the same direction, 
one from each polytope, does not bear any burden on our methods, and neither does the 
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extreme case of two polytopes with identical sets of facet normals. 

The results of experimentation with a broad family of convex polyhedra are reported. 
The relevant programs, source code, data sets, and documentation are available at http : // 
www.cs.tau.ac.il/~efif/CD. A short movie |FH05j that describes some of the concepts 



of the cubical Gaussian-map method can be downloaded from http://acg. cs .tau. a c . 
projects/ internal-projects/ gaussian-map-cubical/Mink3 d. avi[ Another short 
movie [FSHOSbj that describes some of the concepts of the (spherical) Gaussian map method 



among the other can be downloaded from http: //acg. cs .tau. ac . il/projects/ internal- 



[pro j ects/ arr-geodesic-sphere/ mov ie/ aos-xvid. avii 

Both methods are implemented on top of Cgal, and are mainly based on the arrangement 
package of the library (see Chapter [2] and [FWH04| IWFZHOTb] ) . although other parts, such 
as the polyhedral-surface package developed by L. Kettner |Ket99] . are used as well. In 
some cases it is sufficient to build only portions of the boundary of the Minkowski sum 
of two given polytopes to answer collision and proximity queries efficiently. This requires 
locating the corresponding features that contribute to the sought portion of the boundary. 
As both methods we have developed employ two-dimensional arrangements implemented 
on top of the Cgal arrangement package, we harness the ability to answer point-location 
queries efficiently that comes along, to locate corresponding features of two given polytopes. 

The rest of this chapter is organized as follows. The Gaussian maps dual representations 
of polytopes in M.^ are described in Section 13.11 along with some of their properties. In 
Sections 13.21 and 13.31 we show how 3D Minkowski sums can be computed efficiently, after 
the summands are converted to (spherical) Gaussian maps and cubical Gaussian maps, 
respectively. Section 13.41 presents an exact implementation of an efficient collision-detection 
algorithm under translation based on either of the dual representations. In the last section, 
dedicated to experimental results, we highlight 
the performance of our methods on various 
benchmarks. Suggestions for future directions 
are provided in the conclusion chapter in Sec- 
tion [^31 The software access-information along 
with some further design details are provided in 
the Appendix. 



3.1 Gaussian Maps 

The Caussian map G = G{P) of a com- 
pact convex polyhedron P in Euclidean three- 
dimensional space is a set-valued function 
from P to the unit sphere which assigns to 
each point p on the boundary of P the set of out- 
ward unit normals to support planes to P at p. 
Thus, the whole of a facet / of P is mapped un- 
der G to a single point, representing the outward 
unit normal to /. An edge e of P is mapped to 






(c) 



(d) 



Figure 3.1: (a) A tetrahedron, (b) the Gaus- 
sian map of the tetrahedron, (c) a cube, and (d) 
the Gaussian map of the cube. 
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a (geodesic) segment G{e) on S^, whose length is easily seen to be the exterior dihedral angle 
at e. A vertex t> of P is mapped by G to a spherical polygon G{v), whose sides are the image 
under G of edges incident to v, and whose angles are the angles supplementary to the planar 
angles of the facets incident to v; that is, G(ei) and G{e2) meet at angle n — a whenever ci 
and 62 meet at angle a. In other words, G{v) is exactly the "spherical polar" of the link of 
V in P. (The link of a vertex is the intersection of an infinitesimal sphere centered at v with 
P, rescaled, so that the radius is 1.) The above implies that G{P) is combinatorially dual 
to P and an arrangement embedded on the unit sphere [HRS92j . Extending the mapping 
above, by marking each face / = G{v) of the arrangement with its dual vertex v, enables a 
unique inverse Gaussian mapping, denoted by G~^, which maps an extended arrangement 
embedded on the unit sphere back to a polytope boundary. 

An alternative and practical definition follows. A direction in 

can be represented by a point m G Let P be a polytope 
in M^, and let V denote the set of its boundary vertices. For 
a direction u, we define the extremal point in direction u to be 
Xv{u) = aTgrnaXp^v {u, p) , where (■, ■) denotes the inner product. 
The decomposition of §^ into maximal connected regions, so that 
the extremal point is the same for all directions within any region forms the Gaussian map 
of P. For some m G §^ the intersection point of the ray ah emanating from the origin with 
one of the planes hsted below is a central projection of u denoted as Ud, and illustrated on 
the right. The relevant planes are Xd = I, d = 1,2,3, if u lies in the positive respective 
hemisphere, and Xd = —1, d = 1,2,3 otherwise. 

Similar to the Gaussian map, the Cubical Gaussian Map (Cgm) G = G{P) of a polytope 
P in is a set-valued function from P to the six faces of the unit cube whose edges are 
parallel to the major axes and are of length two. A facet / of P is mapped under C to a 
central projection of the outward unit normal to / onto one of the cube faces. Observe that, 
a single edge e of P is mapped to a chain of at most four connected segments that lie in 
four adjacent cube-faces respectively, and a vertex t> of P is mapped to at most five abutting 
convex dual faces that he in five adjacent cube-faces, respectively. The decomposition of the 
unit-cube faces into maximal connected regions, so that the extremal point is the same for 
all directions within any region forms the Cgm of P. Likewise, the inverse Cgm, denoted 
by G~^, maps the six extended arrangement embedded on the six faces of the unit cube to 
the polytope boundary. Figure 13^ shows the Cgm of a tetrahedron. 




3.2 The (Spherical) Gaussian- Map Method 

Armed with the geometry-traits class for geodesic arcs on the sphere (see Section [2.61) . we 
compute Minkowski sums of convex polyhedra, by overlaying their respective Gaussian maps 
represented as arrangements of geodesies on the sphere. Each face / of the Gaussian map 
is extended with the coordinates of its dual vertex v = G~^{f), resulting with a unique 
representation. 
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(a) (b) (c) 

Figure 3.2: (a) A tetrahedron, (b) the Cgm of the tetrahedron, and (c) the Cgm unfolded. Thick 
lines indicate real edges. 

3.2.1 The Representation 

An input model of a polytope is provided as a polyhedral mesh in M^. A polyhedral mesh 
representation consists of an array of boundary vertices and the set of boundary facets, 
where each facet is described by an array of indices into the vertex array. Constructing the 
Gaussian map of a model given in this representation is done indirectly. First, the Cgal 
Polyhedron_3 |Ket99j data-structure that represents the model is constructed. This data 
structure provides quick access to the incidence relations on the polytope features. Then, 
the Gaussian map is constructed from the accessible information stored in the Polyhedron_3 
data-structure. Once the construction of the Gaussian map is complete, the Polyhedron_3 
intermediate representation is discarded. 




Pentagonal Truncated Geodesic Ellipsoid 

Hexecontahedron Icosidodecahedron Sphere level 4 like polyhedron 



Figure 3.3: Gaussian maps of various polytopes. 

The Polyhedron_3 data-structure, like the arrangement DcEL, is based on the imple- 
mentation of an Hds; see Section l273l It consists of extendible vertices, halfedges, and facets 
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and incidence relations on them. It provides methods to traverse all vertices, halfedges, and 
facets, and for local traversals, such as traversing all halfedges incident to a specific vertex. 
It also provides quick access from one halfedge to its twin, and to the incident facet to its left. 
Each vertex and halfedge of the Polyhedron_3 data-structure is extended with a Boolean flag 
that indicates whether the vertex or the halfedge respectively have already been processed 
during the construction of the arrangement that represents the Gaussian map. Each facet 
is extended with the handle of an arrangement vertex. The handle extension of a facet / 
that has already been processed points to the dual vertex v = G{f). The procedure that 
converts an (extended) Polyhedron_3 data-structure P representing a polytope to an (ex- 
tended) Arrangement_on_surf ace_2 data-structure G{P) consists of two steps. First all field 
extensions of all vertices, halfedges, and facets of P are cleared. Then, a recursive function 
provided with an arbitrary vertex w of P as a single parameter is invoked. This function tra- 
verses the halfedges incident to v. When it encounters an unprocessed halfedge e, it obtains 
the normal rii and the vertex-handle extension h{vi) of the facet /i adjacent to the left of e 
and the normal n2 and the vertex-handle extension h{v2) of the facet /2 adjacent to the left 
of the next halfedge in the cyclic chain of halfedges incident to v. Finally, the geodesic short 
arc between ni and n2 is constructed and inserted into the arrangement G{P), as explained 
below, using one of the efficient insertion member-functions of Arrangement_on_surf ace_2; 
see Section r2. 3. 2[ Once the insertion is complete, the halfedge e is marked as processed, v is 
marked as processed once all its incident halfedges are processed. The function recursively 
invokes itself providing an unprocessed vertex adjacent to v, and terminates when no such 
vertex is found. 

Let G indicate the new geodesic arc to be inserted into the arrangement representing the 
Gaussian map. Assume that G is w-monotone with respect to the parameterization defined 
by the geometry-traits class; see Section 12. 6[ That is, it does not intersect the identification 
arc. Let vi, V2, fi, and /2 be the two vertices and two facets as described above. There are 
four cases to handle as follows. 

1. If f 1 and V2 are both null, it implies that this is the first attempt to insert a geodesic 
arc into the arrangement. In this case we call insert_in_f ace_interior(C,/), where 
/ is the single face the DcEL was initialized with; see Section 12. 6[ The handle of the 
two new vertices associated with the endpoints of the newly created geodesic arc G are 
stored in the records of the corresponding facets /i and /2 of P respectively for later 
use. 

2. If Vi is null but V2 is not, we call either insert_f rom_left_vertex(C,f2) or 
insert_f roni_right_vertex(C,V2) depending on whether the existing vertex V2 is to 
the right or to the left of G, and update the vertex-handle field of the corresponding 
facet /i with the new vertex vi. 

3. We handle the analogous case where V2 is null but vi is not similarly. 

4. If both Vi and V2 are not null, we call insert_at_vertices(C,fi ,V2) ■ In this case no 
vertex-handle field needs to be updated. 

If G intersects the identification arc, it is first split at the intersection point into two u- 
monotone arcs Gi and C2. Then, Gi and G2 are inserted according to four cases similar to 
the above. The handling is a bit more intricate. For example, consider the case where vi 
is null but V2 already exists. Assume that Gi reaches the right boundary of the parameter 



50 



Chapter 3. Minkowski Sum Construction 



space and C2 reaches the left boundary (they both meet at an identified point). If V2 is 
associated with the left endpoint of Ci, we first call insert_f rom_lef t_vertex(Ci ,^2) , and 
then insert_f rom_lef t_vertex(C2 where v' is the new vertex associated with the right 
endpoint of Ci introduced while Ci is inserted. Otherwise, V2 must be associated with the 
right endpoint of C2. In this case we first call insert_f rom_right_vertex(C2 ,f2), and then 
insert_f rom_right_vertex(Ci ,f ') , where v' is the new vertex associated with the right 
endpoint of C2 introduced while C2 is inserted. The other three cases are handled similarly. 

We have created a large database of models of polytopes. Table 13.31 lists, for a small 
subset of our polytope collection, the number of features in the arrangement of geodesic arcs 
embedded on the sphere that represents the Gaussian map of each polytope. Recall that the 
number of faces (F) of the Gaussian map is always equal to the number of vertices of the 
polytope. However, the number of vertices (V) of the Gaussian map is either equall to, or 
greater than, the number of facets of the primal representation due to intersections between 
Gaussian-map edges and the identification arc. A similar argument holds for the edges. 
That is, the number of halfedges (HE) of the Gaussian map is either equall to, or greater 
than, twice the number of edges of the primal representation. An edge of the Gaussian 
map that intersects the identification arc must be split at the intersection point into two 
M-monotone geodesic arcs. The table also lists the time in seconds (t) it takes to construct 
the arrangement once the intermediate polyhedron is in place, on a Pentium PC clocked at 
1.7 GHz. 
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(a) 



(b) 



The overlay (see Section 12.4.21 for the exact def- 
inition) of the Gaussian maps of two polytopes 
P and Q respectively identifies all pairs of fea- 
tures of P and Q that have parallel support- 
ing planes, as they occupy the same space on 
the unit sphere, thus, identifying all the pair- 
wise features that contribute to the boundary 
of the Minkowski sum of P and Q. A facet of 
the Minkowski sum is either a facet f oiQ trans- 
lated by a vertex of P supported by a plane par- 
allel to /, or vice versa, or it is a facet parallel 
to two parallel planes supporting an edge of P 
and an edge of Q, respectively. A vertex of the Minkowski sum is the sum of two vertices of 
P and Q respectively supported by parallel planes. 

When the overlay operation progresses, new vertices, edges, and faces of the resulting 
arrangement are created based on features of the two operands. When a new feature is 
created its attributes are updated. There are ten cases that arise and must be handled; see 
Section 12.4.21 for the precise enumeration of the various cases. For example, a new face / is 
induced by the overlap of two faces /i and /2 of the two summands, respectively. The primal 
vertex associated with / is set to be the sum of the primal vertices associated with /i and 
/2, respectively. 



Figure 3.4: (a) The Minkowski sum of a tetra- 
hedron and a cube and (b) the Gaussian map of 
the Minkowski sum. 
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Tetra. © Cube DP © ODP PH © TI E116 © 0E116 



Figure 3.5: Gaussian maps of Minkowski sums. DP — Dioctagonal Pyramid, PH — Pentagonal Hexe- 
contahedron, TI — Truncated Icosidodecahedron, GS4 — Geodesic Sphere level 4, EI16 — Ellipsoid-like 
polyhedron made of 16 latitudes and 32 longitudes. 

Table 13^ lists the number of features (V, HE, F) in the arrangement that represents the 
Gaussian map of the respective Minkowski sums. Table 13751 shows the time in seconds (t) it 
takes to construct the arrangement once the Gaussian maps of the summands are in place. 

3.3 The Cubical Gaussian-Map Method 

While using the Cgm increases the overhead of some operations sixfold, and introduces 
degeneracies that are not present in the case of alternative representations, it simplifies the 
construction and manipulation of the representation, as the partition of each cube face is a 
planar map of segments, a well known concept that has been intensively experimented with 
during recent years. Indeed, all the basic software components that the Cgm layer depends 
on are available in Cgal version 3.3 and higher, while many of the software components 
required by the (spherical) Gaussian map method are expected to appear only in a future 
release of Cgal. The Cgm method, being more mature, exhibits better performance than 
the (spherical) Gaussian map method. One of the reasons for the performance gap is the lack 
of optimized primitives that operate on unnormalized vectors in in case of the (spherical) 
Gaussian map method. Evidently, most of the methods of the geometry-traits class that 
handle geodesic arcs embedded on the sphere project their input w-monotone curves and 
points onto one of the axis-aligned planes every time they are invoked, while all geometric 
objects in case of the Cgm method are projected onto the plane a priori. This creates an 
opportunity for optimization; see Section 16.1.41 In addition, the Cgm data structure, being 
based on components confined to the plane, has a broader recognition, as it can be used in 
restricted environments, e.g., 3D hardware accelerators; see Section l675l 

3.3.1 The Representation 

We use the Cgal Arrangement_2 data-structure to maintain the planar maps. The construc- 
tion of the six planar maps from the polytope features and their incident relations is similar 
to the construction of the arrangement of geodesic arcs that represents the Gaussian map; 
see Section 13.2.11 As in the case of the (spherical) Gaussian map, the Polyhedron_3 inter- 
mediate representation is discarded once the construction of the Cgm is complete. However, 
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while a single edge of P is mapped to at most two u-monotone geodesic arcs in case of the 
(spherical) Gaussian map, it can be mapped to a chain of at most four connected segments 
that lie in four adjacent cube-faces, respectively. In any case the constructions of the Gaus- 
sian maps of both methods respectively amount to the insertion of curves that are pairwise 
disjoint in their interior into the arrangement, an operation that is carried out efficiently, 
especially when one or both endpoints are known. The construction of the Minkowski sum, 
described in Section 13.3.21 amounts to the computation of the six overlays of six pairs of 
planar maps, respectively, an operation well supported by the data structure as well. 

A related dual representation had been considered and discarded before the Cgm rep- 
resentation was chosen. It uses only two planar maps that partition two parallel planes 
respectively instead of six, but each planar map partitions the entire plane. ^ In this 2- map 
representation facets that are near orthogonal to the parallel planes are mapped to points 
that are far away from the origin. The exact representation of such points requires coordi- 
nates with large bit-lengths, which increases significantly the time it takes to perform exact 
arithmetic operations on them. Moreover, facets exactly orthogonal to the parallel planes 
are mapped to points at infinity, and require special handling all together. 

Features that are not in general position, such as two parallel facets facing the same 
direction, one from each pol5d;ope, or worse yet, two identical polytopes, typically require 
special treatment. Still, the handling of many of these problematic cases falls under the "gen- 
eral" case, and becomes transparent to the Gaussian-map layer (either cubical or spherical). 
Consider for example the case of two neighboring facets in one polytope that have parallel 
neighboring facets in the other. This translates to overlapping segments in case of the Cgm 
method and overlapping geodesic arcs in case of the (spherical) Gaussian-map method, one 
from each Gaussian map of the two polytopes,^ that appear during the Minkowski sum com- 
putation. The algorithm that computes it is oblivious to this condition, as the underlying 
arrangement data structure, either embedded in the plane or on the sphere, is perfectly ca- 
pable of handling overlapping curves. However, as mentioned above, other degeneracies do 
emerge, and are handled successfully. One example, in case of the Cgm, is a facet / mapped 
to a point that lies on an edge of the unit cube, or even worse, coincides with one of the eight 
corners of the cube. Figure I3.7( a.b.c) depicts an extreme degenerate case of an octahedron 
oriented in such a way that its eight facets are mapped to the eight vertices of the unit cube, 
respectively. The (spherical) Gaussian map method is not free of degenerate conditions. For 
example, a facet mapped to a point that lies on the identification arc, or worse yet, coincides 
with one of the two poles. 

The dual representation is extended further, in order to handle all these degeneracies and 
perform all the necessary operations as efficiently as possible. Each planar map is initialized 
with four edges and four vertices that define the unit square — the parallel-axis square 
circumscribing the unit circle. During construction, some of these edges or portions of them 
along with some of these vertices may turn into real elements of the Cgm. The introduction 



^Each planar map that corresponds to one of the six unit-cube faces in the Cgm representation also 
partitions the entire plane, but only the [—1,-1] x [1,1] square is relevant. The unbounded face, which 
comprises all the rest, is irrelevant. 

^Other conditions translate to overlapping segments in case of the Cgm or geodesic arcs in case of the 
(spherical) Gaussian map method as well. 
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of these artificial elements not only expedites the traversals below, but is also necessary for 
handling degenerate cases, such as an empty cube face that appears in the representation 
of the tetrahedron and depicted in Figure [221(c). The global data consists of the six planar 
maps and 24 references to the planar vertices that coincide with the unit-cube corners. 

The exact mapping from a facet nor- 
mal in the 3D coordinate-system to a pair 
that consists of a planar map and a planar 
point in the 2D coordinate-system is de- 
fined precisely through the indexing and 
ordering system, illustrated in Figure 13.61 
Now before your eyes cross permanently, 
we advise you to keep reading the next 
few lines, as they reveal the meaning of 
some of the enigmatic numbers that ap- 
pear in the figure. The six planar maps 
are given unique ids from through 5. Ids 
0, 1, and 2 are associated with the planes 
X = —1, y = —1, and z = —1, respectively, 
and ids 3, 4, and 5 are associated with the 
planes x = 1, y = 1, and z = 1, respec- 
tively. The major axes in the 2D Cartesian 
coordinate-system of each planar map are 
determined by the 3D coordinate-system. 
The four corner vertices of each planar 
map are also given unique ids from through 3 in lexicographic order in their respective 2D 
coordinate-system, see Table \37L\ columns titled Underlying Plane and 2D Axes. 

Table 3.1: The coordinate systems and the cyclic chains of corner vertices. PM stands for Planar 
Map, and Cr stands for Corner. 



Underlying 


2D Axes 


Corner 




Plane 


(0,0) 


1 (0,1) 


2 (1,0) 


3 (1,1) 


Id 


Eq 


X 


Y 


PM 


Cr 


PM 


Cr 


PM 


Cr 


PM 


Cr 





X = —1 


z 


Y 


1 





2 


2 


5 





4 


2 


1 


y = -l 


X 


Z 


2 








2 


3 





5 


2 


2 


z = -l 


Y 


X 








1 


2 


4 





3 


2 


3 


X = 1 


Y 


Z 


2 


1 


1 


3 


4 


1 


5 


3 


4 


y = l 


Z 


X 





1 


2 


3 


5 


1 


3 


3 


5 


z = 1 


X 


Y 


1 


1 





3 


3 


1 


4 


3 



Each feature type of the DcEL used to maintain the incidence relations of the vertices, 
halfedges, and faces of the Arrangement_2 data structure (see Section 12.31) is extended to 
hold additional attributes. Some of the attributes are introduced only in order to expedite 
the computation of certain operations, but most of them are necessary to handle degenerate 
cases such as a planar vertex lying on the unit-square boundary. Each planar-map vertex 



A' 




Y 



Figure 3.6: The data structure. Large-font numbers 
indicate plane ids. Small-font numbers indicate corner 
ids. X and Y axes in different 2D coordinate systems 
are rendered in different colors. 
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V is extended with (i) the coefficients of the plane containing the polygonal facet C~^{v) 
(see Section 13.11 for the definition of C and C~^), (ii) the location of the vertex — an 
enumeration indicating whether the vertex coincides with a cube corner, or lies on a cube 
edge, or contained in a cube face, (iii) a Boolean flag indicating whether it is non-artificial 
(there exists a facet that maps to it), and (iv) a pointer to a vertex of a planar map associated 
with an adjacent cube-face that represents the same central projection for vertices that 
coincide with a cube corner or lie on a cube edge. Each planar-map halfedge e is extended 
with a Boolean flag indicating whether it is non-artificial (there exists a polytope edge that 
maps to it). Each planar- map face / is extended with the polytope vertex that maps to it 
v = C-\f). 




(d) (e) (f) 




(g) (h) (i) 

Figure 3.7: (a) An octahedron, (d) a dioctagonal pyramid, (g) an ellipsoid-like polyhedron level 16, 
(b,e,h) the Cgm of the respective polytope, and (c,f,i) the Cgm unfolded. 

Each vertex that coincides with a unit-cube corner or lies on a unit-cube edge contains 
a pointer to a vertex of a planar map associated with an adjacent cube face that represents 
the same central projection. Vertices that lie on a unit-cube edge (but do not coincide with 
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1,0 



2,0 



unit-cube corners) come in pairs. Two vertices that form such a pair he on 
the unit-square boundary of planar maps associated with adjacent cube 
faces, and they point to each other. Vertices that coincide with unit- 
cube corners come in triplets and form cyclic chains ordered clockwise 
around the respective vertices. The diagram on the right specifies one of 
the eight cyclic chains. The left and right indices of each pair specify a planar-map id and 
a corner id, respectively. All specific connections are listed in Table 13.11 Recall that all 
maps are facing outwards. As a convention, edges incident to a vertex are ordered clockwise 
around the vertex, and edges that form the boundary of a face are ordered counterclockwise. 
The Polyhedron_3 and Arrangement_2 data structures for example both use a DCEL data 
structure that follows the convention above. 

We provide a fast clockwise traversal of the faces incident to any given 
vertex v. Clockwise traversals around internal vertices are immediately 
available by the DcEL. Clockwise traversals around boundary vertices 
are enabled by the cyclic chains above. This traversal is used to calculate 
the normal to the (primary) polytope-facet / = C~^{v) and to render 
the facet. Fortunately, rendering systems are capable of handling a 
sequence of vertices that define a polygon in clockwise order as well, an 
order opposite to the conventional ordering above. 

The data structure also supports a fast traversal over the planar-map 
halfedges that form each one of the four unit-square edges. This traversal 
is used during construction to quickly locate a vertex that coincides with a 
cube corner or lies on a cube edge. It is also used to update the cyclic chains 
of pointers mentioned above; see Section 13.3.21 

We maintain a flag that indicates whether a planar vertex coincides with a cube corner, 
a cube edge, or a cube face. At first glance this looks redundant. After all, this information 
could be derived by comparing the x and y coordinates to —1 and -|-1. However, it has 
a good reason as explained next. Using exact number-types often leads to representations 
of the geometric objects with large bit-lengths. Even though we use various techniques to 
prevent the length from growing exponentially [FWH04] , we cannot prevent the length from 
growing at all. Even the computation of a single intersection requires a few multiplications 
and additions. Cached information computed once and stored at the features of the planar 
map avoids unnecessary processing of potentially long representations. 





The table to the right shows the number of ver- 
tices (V), halfedges (HE), and faces (F) of the six 
planar maps that comprise the Cgm of the dioctago- 
nal pyramid shown in Figure 13771 (d,e,f). The number 
of faces of each planar map include the unbounded 
face. Table 13.31 shows the number of features in the 
primal and dual representations of a small subset of 
our polytopes collection, on which we report the re- 
sults of experiments below. The number of planar 
features is the total number of features of the six pla- 
nar maps. 



Table 3.2: The number of features of 
the six planar maps of the Cgm of the 
dioctagonal pyramid. 



Planar map 


V 


HE 


F 


0, (x = -1) 


12 


32 


6 


1, (2/ = -l) 


28 


80 


14 


2, (^ = -1) 


12 


32 


6 


3, (x = 1) 


12 


32 


6 


4, {y = 1) 


21 


72 


17 


5, {z = 1) 


12 


32 


6 


Total 


97 


280 


55 
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3.3.2 Exact Minkowski Sums 

A similar argument regarding the representation of Minkowski sums using Gaussian maps 
mentioned in Section 13.21 holds for the cubical Gaussian maps with the unit cube replacing 
the unit sphere. More precisely, a single map that subdivides the unit sphere is replaced 
by six planar maps, and the computation of a single overlay is replaced by the computation 
of six overlays of corresponding pairs of planar maps. Recall that each (primal) vertex is 
associated with a planar-map face, and is the sum of two vertices associated with the two 
overlapping faces of the two Cgm's of the two input polytopes, respectively. 

Each planar map in a Cgm is a convex subdivision. Finke and Hinrichs |FH95j describe 
how to compute the overlay of such special subdivisions optimally in linear time. However, 
a preliminary investigation shows that a large constant governs the linear complexity, which 
renders this choice less attractive. Instead, we resort to a sweep-line based algorithm that 
exhibits good practical performance, and incurs a mere logarithmic factor over the optimal 
computing time. In particular we use the overlay operation supported by the Arrangement 
package. It requires the provision of a complementary component that is responsible for up- 
dating the attributes of the DcEL features of the resulting six planar maps; see Section [2.4.21 

The overlay operates on two instances of Arrangement_2. In the description below Vi, 
Ci, and /i denote a vertex, a halfedge, and a face of the first operand respectively, and V2, 
62, and /2 denote the same feature types of the second operand, respectively. When the 
overlay operation progresses, new vertices, halfedges, and faces of the resulting planar map 
are created based on features of the two operands. Exactly ten cases described below arise 
and must be handled. When a new feature is created its attributes are updated. The updates 
performed in all cases except for case (1) are simple and require constant time. 

1. A new vertex v is induced by coinciding vertices Vi and V2- 

The location of the vertex v is set to be the same as the location of the vertex vi 
(the locations of V2 and vi must be identical). The induced vertex is not artificial if 
and only if (i) at least one of the vertices vi or V2 is not artificial, or (ii) the vertex 
lies on a cube edge or coincides with a cube corner, and both vertices Vi and V2 have 



Table 3.3: Complexities of the primal and dual representations. DP — Dioctagonal Pyramid, PH — 
Pentagonal Hexecontahedron, TI — Truncated Icosidodecahedron, GS4 — Geodesic Sphere level 4, EI16 
— Ellipsoid-like polyhedron made of 16 latitudes and 32 longitudes, t - time consumption in seconds. 



Object Type 


Primal 


SGM 


CGM 


V 


E 


F 


V 


HE 


F 


t 


V 


HE 


F 


t 


Tetrahedron 


4 


6 


4 


4 


12 


4 


0.01 


42 


102 


21 


0.01 


Octahedron 


6 


12 


8 


10 


28 


6 


0.01 


24 


48 


12 


0.01 


Icosahedron 


12 


30 


20 


21 


62 


12 


0.01 


72 


192 


36 


0.01 


DP 


17 


32 


17 


25 


80 


17 


0.01 


97 


280 


55 


0.01 


PH 


60 


150 


92 


101 


318 


60 


0.03 


200 


600 


112 


0.02 


TI 


120 


180 


62 


77 


390 


120 


0.05 


230 


840 


202 


0.03 


GS4 


252 


750 


500 


506 


1512 


252 


0.08 


708 


2124 


366 


0.07 


E116 


482 


992 


512 


528 


2016 


482 


0.11 


776 


2752 


612 


0.06 
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non-artificial incident halfedges that do not overlap. 

2. A new vertex v is induced by a vertex vi that lies on an edge 62- 

The location of the vertex v is set to be the same as the location of the vertex Vi. v is 
not artificial if and only if vi is not artificial or 62 is not artificial. 

3. An analogous case of a vertex V2 that lies on an edge ei. 

4. A new vertex v is induced by a vertex Vi that is contained in a face /2. 

The attributes of the vertex v are set to be the same as the attributes of the vertex vi . 

5. An analogous case of a vertex V2 contained in a face /i. 

6. A new vertex v is induced by the intersection of two edges Ci and 62- 

The vertex v cannot lie on a cube edge and cannot coincide with a cube corner. Thus, 
it is necessarily not artificial. 

7. A new edge e is induced by the overlap of two edges ei and 62- 
The edge e is not artificial if at least one of Ci or 62 is not artificial. 

8. A new edge e is induced by an edge ei that is contained in a face f2- 
The edge e is not artificial if ei is not artificial. 

9. An analogous case of an edge 62 contained in a face /i. 

10. A new face / is induced by the overlap of two faces fi and f2- 

The primal vertex associated with / is set to be the sum of the primal vertices associ- 
ated with fi and /2, respectively. 

After the six map overlays are computed, some maintenance operations must be per- 
formed to obtain a valid Cgm representation. As mentioned above, the global data consists 
of the six planar maps and 24 references to vertices that coincide with the unit-cube cor- 
ners. For each planar map we traverse its vertices, obtain the four vertices that coincide 
with the unit-cube corners, and initialize the global data. We also update the cyclic chains 
of pointers to vertices that represent identical central projections. To this end, we exploit 
the fast traversal over the halfedges that coincide with the unit-cube edges mentioned in 
Section 13.3. 11 

The complexity of a single overlay operation is 0{klogn), where n is the total number 
of vertices in the input planar maps, and k is the number of vertices in the resulting planar 
map. The total number of vertices in all the six planar maps in a Cgm that represents a 
polytope P is of the same order as the number of facets in the primary polytope P. Thus, 
the complexity of the entire overlay operation is 0{F\og{Fi + F2)), where Fi and F2 are the 
number of facets in the input polytopes respectively, and F is the number of facets in the 
Minkowski sum. 



3.4 Exact Collision Detection 

Computing the separation distance between two polytopes with m and n features respec- 
tively can be done in O(logmlogn) time, after an investment of at most linear time in 
preprocessing [DK90j . Many practical algorithms that exploit spatial and temporal coher- 
ence between successive queries have been developed, some of which became classic, such as 
the GJK algorithm jGJK88j and its improvement |Cam97j . and the LC algorithm |LC91] and 
its optimized variations [ELOOl IGHZ99i IMir98j . Several general-purpose software libraries 
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that offer practical solutions are available today, such as the SOLID library [21] based on 
the improved GJK algorithm, the SWIFT library [26] based on an advanced version of the 
LC algorithm, the QuickCD library [21] , and more. For an extensive review of methods and 
libraries see the survey [LM04j . 

Given two polytopes P and Q, detecting collision between them and computing their rel- 
ative placement can be conveniently done in the configuration space, where their Minkowski 
sum M = P(B{—Q) resides. These problems can be solved in many ways, and not all require 
the explicit representation of the Minkowski sum M. However, having it available is attrac- 
tive, especially when the polytopes are restricted to translations only, as the combinatorial 
structure of the Minkowski sum M is invariant to translations of P or Q. The algorithms 
described below are based on the following well known observations (see Chapter [1] for defi- 
nitions): 

p«nQ'"^0^w-nGM = P© i-Q) , 
7r(P",g'") =min{||t|| \ {w - u + t) e M,t e M.^} , 
5,(P", Q"") = mi{a \{w -u + af) ^ M,a eR} . 

Given two polytopes P and Q in either (spherical) Gaussian map or Cgm representation 
respectively, we reflect Q through the origin to obtain —Q, compute the Minkowski sum 
M = P © {—Q), and retain it in the respective Gaussian-map representation G{M). Then, 
each time P or Q or both translate by two vectors u and w in respectively, we apply a 
procedure that determines whether the query point s = w—u is inside, on the boundary of, or 
outside M. In addition to an enumeration of one of the three conditions above, the procedure 
returns a witness of the respective relative placement. Let r be a ray emanating from an 
internal point c G M and going through s. If the (spherical) Gaussian map representation 
is used, the witness data is the vertex v = G{f) — a mapping of a facet f of M embedded 
on the sphere and stabbed by the ray r. If the Cgm representation is used, the witness 
data is a pair that consists of a vertex v = G{f) — a mapping of a facet f of M embedded 
in a unit cube face, and the planar map V containing v. This information is used as a 
hint in consecutive invocations. The internal point c could be the average of all vertices 
of M computed once and retained along M, or just the midpoint of two vertices that have 
supporting planes with opposite normals easily extracted from either map representation. 
Once / is obtained, determining whether P" and collide is trivial, according to the first 
formula (of the three) above. The query point s is contained in the open half-space defined 
by the supporting plane to / if and only if s is outside of M, this occurs if and only if P" 
does not collide with Q^. 

The collision-detection procedure applies a local walk 
on the respective Gaussian map faces. It starts with some 
vertex vq, and then performs a loop moving from the cur- 
rent vertex to a neighboring vertex, until it reaches the 
final vertex. If the Cgm representation is used, the pro- 
cedure may jump from a planar map associated with one 
cube-face to a different one associated with an adjacent 
cube-face. The first time the procedure is invoked, vq is 
chosen to be a vertex that lies on the central projection of Figure 3.8: Simulation of motion. 
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the normal directed in the same direction as the ray r. In consecutive calls, vq is chosen 
to be the final vertex of the previous call exploiting spatial and temporal coherence. The 
figure above is a snapshot of a simulation program that detects collision between a static 
obstacle and a moving robot, and draws the obstacle and the trail of the robot; instructions 
to obtain, install, and execute the program appear in the Appendix. The Minkowski sum is 
recomputed only when the robot is rotated, which occurs every other frame. The program 
has the distinctive feature of being able to identify the case where the robot grazes the ob- 
stacle, but does not penetrate it, since it produces exact results. The computation takes 
just a fraction of a second on a Pentium PC clocked at 1.7 GHz using either representation. 
Similar procedures that compute the directional penetration-depth and minimum distance 
are available as well. 
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(a) (b) (c) 

Figure 3.9: (a) The Minkowski sum (of two polytopes) the complexity of which is maximal, (b) the 
Cgm of the Minkowski sum, and (c) the Cgm unfolded. Red lines are graphs of edges that originate 
from one polytope and blue lines are graphs of edges that originate from the other. 

In Chapter H] we show that the exact maximum number of facets of Minkowski sums of 
two polytopes is 4mn — 9m — 9n -|- 26, where m and n are the number of facets of the two 
summands, respectively. This bound is tight |FHW07] . The example depicted in Figure 13.91 
shows a Minkowski sum that reaches this maximum complexity. It is the sum of two identical 
polytopes, each containing n faces (n = 11 in Figure 13.91) . but one is rotated about the 
vertical axis approximately^ 90° relative to the other. The polytopes are specifically shaped 
to ensure that the number of intersections between dual edges, which are the mappings of 
the polytope edges, is maximal. A careful counting reveals that the number of vertices in the 
dual representation excluding the artificial vertices reaches 4-11-11 — 9-11 — 9-11 + 26 = 312, 
which is the number of facets of the Minkowski sum. 

Not every pair of polytopes yields a Minkowski sum proportional to mn. As a matter 
of fact, it can be as low as n in the extremely-degenerate case of two identical polytopes 
variant under scaling. Even if no degeneracies exist, the complexity can be proportional to 



■^The results of all rotations are approximate, as we have not yet dealt with exact rotation. One of our 
future goals is the handling of exact rotations. 
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only m + n, as in the case of two geodesic spheres^ level / = 2 slightly rotated with respect 
to each other, depicted in Figure I3.10[ Naturally, an algorithm that accounts for all pairs 
of vertices, one from each polytope, is rendered inferior compared to an output-sensitive 
algorithm like ours in such cases, as we demonstrate in the next section. 




(a) (b) (c) 

Figure 3.10: (a) The Minkowski sum of two geodesic spheres level 2 slightly rotated with respect to 
each other, (b) the Cgm of the Minkowski sum, and (c) the Cgm unfolded. 

3.6 Experimental Results 




(a) (b) (c) 

Figure 3.11: (a) The Minkowski sum of two approximately orthogonal squashed dioctagonal pyramids, 
(b) the Cgm, and (c) the Cgm unfolded. 

As mentioned above, the Minkowski sum of two polytopes is the convex hull of the pair- 
wise sum of the vertices of the two polytopes. We have implemented this straightforward 
method using the Cgal convex_hull_3 function, which uses the Polyhedron_3 data struc- 
ture to represent the resulting polytope, and used it to verify the correctness of our two 
methods. We compared the time it took to compute exact Minkowski sums using our two 
methods, a third method implemented by Hachenberger based on Nef polyhedra embed- 
ded on the sphere |HKM07] . a fourth method implemented by Weibel [28], based on an 
output-sensitive algorithm designed by Fukuda [Fuk04j . and the naive convex-hull method. 

"^An icosahedron, every triangle of which is divided into 4' triangles using class I aperture 4 partition 
method, whose vertices are elevated to the circumscribing sphere (S WK03j . 
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The Nef-based method is not speciahzed for Minkowski sums. It can compute the overlay 
of two arbitrary Nef polyhedra embedded on the sphere, which can have open and closed 
boundaries, facets with holes, and lower dimensional features. The overlay is computed by 
two separate hemisphere-sweeps. 

Fukuda's algorithm relies on linear programming. Its complexity is 0(5LP(3, 5)V), where 
5 = 5i + ^2 is the sum of the maximal degrees of vertices, 5i and 82, in the two input polytopes 
respectively, V is the number of vertices of the resulting Minkowski sum, and LP{d, m) is 
the time required to solve a linear programming in d variables and m inequalities. Note 
that Fukuda's algorithm is more general, as it can be used to compute the Minkowski sum 
of polytopes in an arbitrary dimension d, and as far as we know, it has not been optimized 
specifically for d = 3. 

Table 3.4: Complexities of primal and dual Minkowski-sum representations. DP — Dioctagonal Pyra- 
mid, ODP — Dioctagonal Pyramid orthogonal to DP, PH — Pentagonal Hexecontahedron, TI — 
Truncated Icosidodecahedron, GS4 — Geodesic Sphere level 4, RGS4 — Rotated Geodesic Sphere 
level 4, EI16 — Ellipsoid-like polyhedron made of 16 latitudes and 32 longitudes, 0EI16 — Ellipsoid-like 
polyhedron made of 16 latitudes and 32 longitudes orthogonal to EI16. 
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Table 3.5: Time consumption (in seconds) of the Minkowski-sum computation. CH — the Convex 
Hull method, SGM — the (spherical) Gaussian map based method, CGM — the cubical Gaussian-map 
based method, NGM — the Nef based method, Fuk — Fukuda's Linear-Programming based algorithm, 
^^p"^ — the ratio between the product of the number of input facets and the number of output facets. 
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E116 
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The results listed in Table 13.51 produced by experiments conducted on a Pentium PG 
clocked at 1.7 GHz, show that our methods are much more efficient in all cases, and the 
Cgm method in particular is more than fifty times faster than the convex-hull method in 
one case. The number of models used as summands in the listed experiments is just a small 
fraction of the total number of models in our collection, which contains hundreds of models 
of polytopes; see Section IA.2I for instructions how to download the corresponding files. The 
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listed experiments are just a small sample of all the experiments we have conducted. The 
last column of the table indicates the ratio ^^y^, where Fi and F2 are the number of facets 
of the input polytopes respectively, and F is the number of facets of the Minkowski sum. As 
this ratio increases, the relative performance of the output-sensitive algorithms compared to 
the convex-hull method improves as expected. Figure 13.121 illustrates some of the resulting 
Minkowski sums listed in Table 1331 

The SGM, CGM, NGM, and CH methods are all based on Cgal. As the implementa- 
tion of these methods and of Cgal is generic, it is possible to instantiate specific components, 
such as the number type, the geometric kernel, and the segment-handling traits class with 
the model of the respective concept that achieves the best result. The code of the programs 
used to obtain the results listed in Table 13.51 are instantiated with the CGAL : : Gmpq exact 
rational number-type, the Simple_cartesian geometric kernel, and the Lazy_kernel kernel 
adaptor, which performs lazy exact computations [PF06j . The computation is delayed until 
a point where the approximation with interval arithmetic is not precise enough to perform 
safe comparisons. In other words, these programs need only to compute to sufficient preci- 
sion to evaluate predicates correctly, exploiting a significant relaxation of the naive concept 
of numerical exactness. 

We experimented with two different exact number types: One provided by Leda 4.4.1, 
namely leda: :rational, and another based on GMP 4.1.2, namely CGAL: :Gmpq. The 
former does not normalize the rational numbers automatically. Therefore, we had to initiate 
normalization operations to contain their bit-length growth. In case of the Cgm method, 
for example, we chose to do it right after the central projections of the facet-normals are 
calculated, and before the chains of segments, which are the mapping of facet-edges, are 
inserted into the planar maps. Our experience shows that indiscriminate normalization 
considerably slows down the planar-map construction, and the choice of number type may 
have a drastic impact on the performance of the code overall. 
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(d) (e) (f) 




(g) (h) (i) 



Figure 3.12: (a) The Minkowski sum of two approximately orthogonal dioctagonal pyramids, (d) 
the Minkowski sum of a Pentagonal Hexecontahedron and a Truncated Icosidodecahedron, (g) the 
Minkowski sum of two approximately orthogonal ellipsoid-like polyhedra, (b,e,h) the Cgm of the re- 
spective polytope, and (c,f,i) the Cgm unfolded. 
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/ was working on the proof 
of one of my poems all the 
morning, and took out a 
comma. In the afternoon I 
put it back again. 

Oscar Wilde 




Exact Complexity of Minkowski Sums 



We present a tight bound on the exact maximum complexity of Minkowski sums of polytopes 
in M^. In particular, we prove that the maximum number of facets of the Minkowski sum of k 
polytopes with mi, m2, . . . , m^ facets respectively is bounded from above by J2i<i<j<ki'^^i~ 
5)(2mj — 5) + X]i<j<A: + (2)- Given k positive integers mi, m2, . . . , m^, we describe how to 
construct two polytopes with corresponding number of facets, such that the number of facets 
of their Minkowski sum is exactly X]i<j<j<A:(2"^i ~ 5){2mj — 5) + X]i<j<A:"^i + (2)- When 
k = 2, for example, the expression above reduces to 4mim2 — 9mi — 9m2 + 26. Figure H]T] 
illustrates some polytopes that when rotated properly and used as summand, the number of 
facets of the resulting Minkowski sums is maximal. 

Various methods to compute the Minkowski sum of two polyhedra in have been 
proposed; see Section 11.21 for details about these methods and about the combinatorial 
complexity of the sum. Recall, that (i) a common approach for computing Minkowski sums 







(b) 



(c) 



(d) 



Figure 4.1: Gaussian maps of summands of Minkowski sums with maximal number of facets. The 
polytopes represented by the Gaussian maps (a), (b), (c), and (d) consist of 4, 5, 11, and 101 facets, 
respectively. 
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of non-convex polyhedra decomposes each polyhedron into convex pieces, and computes 
pairwise Minkowski sums of the convex pieces, and (ii) all the efficient methods are output 
sensitive. Thus, the exact maximum complexity of the Minkowski sum structure has practical 
implications. 

One method to compute the Minkowski sum of two polytopes is to compute the convex 
hull of the pairwise sum of the vertices of the two polytopes. While being simple and easy 
to implement, the time complexity of this method is Q{mn) regardless of the size of the 
resulting sum, which can be as low as (m + n) (counting facets) for degenerate cases. ^ In 
Chapter [3] we describe several complete implementations of output-sensitive methods for 
computing exact Minkowski sums (beyond the naive method mentioned above), including 
two methods that we introduce. These methods exploit efficient innovative techniques in 
the area of exact geometric computing to minimize the time it takes to ensure exact results. 
However, even with the use of these techniques, the amortized time of a single arithmetic 
operation is larger than the time it takes to carry out a single arithmetic operation on 
native number types, such as floating point. Thus, the constant that scales the dominant 
element in the expression of the time complexity of these algorithms increases, which makes 
the question this chapter attempts to answer, "What is the exact maximum complexity of 
Minkowski sums of polytopes in R^?", even more relevant. 

Gritzmann and Sturmfels [GS93j formulated an upper bound on the number of features 
ff of any given dimension i of the Minkowski sum of many polytopes in d dimensions: 

d-i-l 

ff{Pi © A © • • • © ^fe) < 2 (^) J2 (^X^) for 2 = 0, 1, ... c/ - 1, where j denotes the number 

h=0 

of non-parallel edges of Pi, P2, Pk- According to this expression, the number of facets /| 
of the Minkowski sum of two polytopes in is bounded from above by j(j — 1). Fukuda and 
Weibel [FW07J obtained upper bounds on the number of edges and facets of the Minkowski 
sum of two polytopes in M.^ in terms of the number of vertices of the summands: f2{Pi®P2) < 
/q (Pi)/q (P2) + fo{Pi) + /o (-^2) — 6. They also studied the properties of Minkowski sums of 
perfectly centered polytopes and their polars, and provided a tight bound on the number of 
vertices of the sum of polytopes in any given dimension. 

The main result presented in this chapter concerning two polytopes follows. 

Theorem 4.1. Let Pi, P2, . . . ,Pk be a set ofk polytopes in M^, such that the number of facets 
of Pi is rrii for i = 1,2, . . . , k. The number of facets of the Minkowski sum Pi © P2 © • • • © -Pfc 

cannot exceed J2i<i<:j<ki'^''^i~^){'^^j~^)+J2'i=i^i~^ {2) ■ This bound is tight. Namely, given 
k integers mi, m2, ■ ■ ■ , rn^, such that rrii > A for i = 1,2, . . . , k, it is possible to construct k 
polytopes in with corresponding number of facets, such that the number of facets of their 
Minkowski sum is exactly the expression above. 

The rest of this chapter is organized as follows. The upper bound on the number of facets 
of Minkowski sums for the special case of two polytopes in is derived in Section 14.11 In 
Section 14.21 we describe how to construct two polytopes, the number of facets of which 
is given, such that the number of facets of their Minkowski sum is identical to the bound 
derived in Section I^TTl The bounds for the general case of k polytopes is proved in Section 



^It can be as low as m(= n) in the extremely-degenerate case of two similar polytopes with parallel facets. 
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Information about the polyhedra models and the interactive program that computes their 
Minkowski sums and visuahzes them, used to verify the results and generate the figures in 
this chapter are provided in the Appendix. 

4.1 The Upper Bound for k = 2 

The overlay (see Section [2.4.21 for the exact definition) of the Gaussian maps (see Section ISTTl 
for the exact definition) of two polytopes P and Q respectively is the Gaussian map of the 
Minkowski sum of P and Q; see Section [3.2.21 for a detailed explanation. 

The number of facets of the Minkowski sum M of two polytopes P and Q with m and n 
facets respectively is equal to the number of vertices of the Gaussian map G{M) of M. A 
vertex in G{M) is either due to a vertex in the Gaussian map of P, or due to a vertex in 
the Gaussian map of Q, or due to an intersection of exactly two edges, one of the Gaussian 
map of P and the other of the Gaussian map of Q. Thus, the number of facets of M cannot 
exceed m + n + g{M), where g{M) is the number of intersections of edges of G{P) with edges 
of G(Q) in G(M).2 

Observation 4.2. The maximum exact number of edges in a Gaussian map G{P) of a 
polytope P with m facets is 3m — 6. The maximum exact number of faces is 2m — 4. Both 
maxima occur at the same Gaussian maps. 

The above can be obtained by a simple application of Euler's formula for planar graphs 
to the Gaussian map G{P). It can be used to trivially bound the exact maximum number of 
facets of the Minkowski sum of two polytopes defined as /(m, n) = max{/(P © Q) \ f{P) = 
'^ifiQ) = ''T'}: where f{P) is the number of facets of a polytope P. First, we can use 
the bound on the number of edges to obtain: f{m,n) < m + n + (3m — 6) ■ (3n — 6) = 
9mn — 17m — 17n + 36. Better yet, we can plug the bound on the number of dual faces, 
which is the number of primal vertices, in the expression introduced by Fukuda and Weibel, 
see above, to obtain: f{m,n) < {2m—A)-{2n~A) + {2m—4) + {2n—A)—6 = Amn—6m—6n+2. 
Still, we can improve the bound even further, but first we need to bound the number of faces 
in G{M). 

Lemma 4.3. Let Gi and G2 be two Gaussian maps of convex polytopes, and let G be their 
overlay. Let fi, f2, and f denote the number of faces of Gi, G2, and G, respectively. Then, 

/</l-/2. 

Each face in the overlay is an intersection of a face of each map. Since these faces are 
spherically convex (and smaller than hemispheres), the intersection is also spherically convex 
(and thus connected). This lemma is similar to the one where convex planar maps replace 
the Gaussian maps. Nevertheless, we provide a formal proof directly applied to the spherical 
case. 

Proof. We label each face in G by a pair of indices of the originating overlaid faces in Gi 
and G2 respectively, and argue that no two faces in G can have the same label. Assume to 



^The number of facets is strictly equal to the given expression, only when no degeneracies occur. 
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the contrary that there exist two faces ha and hf, in G that have the same label, say {i,j)- 
That is, the faces h{ in Gi and h^ in G2 induce the two distinct faces ha and hb- Pick two 
points a E ha and h E h^. There must be a geodesic segment between a and h that is entirely 
contained in h\ and also in /i2; both maps are spherically convex. This implies that none 
of the edges in Gi and G2 split this geodesic segment, contradicting the fact that they reside 
in two different faces of G. □ 

We are ready to tackle the upper bound of Theorem 14.11 for the special case k = 2, that 
is, prove that the number of facets of the Minkowski sum P © Q of two polytopes P and Q 
with m and n facets respectively cannot exceed Amn — 9m — 9n + 26; see Page [661 

Proof. Let vi, ei, /i and V2, 62, /2 denote the number of vertices, edges, and faces of G{P) and 
G{Q), respectively. The number of vertices, edges, and faces of G{M) is denoted as v, e, and 
/, respectively. Assume that P and Q are two polytopes, such that the number of facets of 
their Minkowski sum is maximal. Recall that the number of facets of a polytope is equal to 
the number of vertices of its Gaussian map. Thus, we have vi = m, V2 = n, and v = f{m, n). 
First, we need to show that vertices of G{P), vertices of G{Q), and intersections between 
edges of G{P) and edges of G{Q) do not coincide. Assume to the contrary that some do. 
Then, one of the polytopes P or Q or both can be slightly rotated to escape this degeneracy, 
but this would increase the number of vertices v = f{m,n), contradicting the fact that 
/(m, n) is maximal. Therefore, the number of vertices v is exactly equal to vi + V2 + Vx, 
where Vx denotes the number of intersections of edges of G{P) and edges of G{Q) in G{M). 
Counting the degrees of all vertices in G{M) implies that 2ei + 2e2 + 4:Vx = 2e. Using Euler's 
formula, we get ei + 62 + 2vx = f + Vi + V2 + Vx — 2. Applying Lemma 14. 3[ we can bound 
Vx from above Vx < /1/2 + fi + ^2 — 2 — ei — 62- 

Observation 14. 21 sets an upper bound on the number of edges ei. Thus, ei can be expressed 
in terms of £1, a non- negative integer, as follows: ei = 3vi —6 — ii. Applying Euler's formula, 
the number of facets can be expressed in terms of ii as well: /i = ei + 2 — f 1 = 2vi — 4 — £1. 
Similarly, we have 62 = 3v2 — Q — £2 and /2 = 2f 2 — 4 — ^2 for some non- negative integer £2- 

Vx < {2vi - 4 - £i){2v2 -A-i2) + vi+V2-2- {3vi - 6 - £1) - {3v2 - Q - £2) 

< AV1V2 - lOvi - Wv2 + 26 + h{£i,£2) , (4.1) 

where 4) = ^1^2 + 5£i + 5^2 - 2vi£2 - 2v2£i. 

G{P) consists of a single connected component. Therefore, the number of edges ei must 
be at least vi — 1. This is used to obtain an upper bound on £1 as follows: vi — 1 < ei = 
3vi —6 — ^1, which implies £1 < 2vi — 5, and similarly £2 < 2v2 — 5. Thus, we have: 

h{£i,£2) = £i£2 + 5£i + M2 - 2vi£2 - 2v2£i 

= ^i(|-(2^2-5)) + £2(|-(2t;i-5))<0. 

From Equation (14. ip we get that Vx < 4fif2 — lOfi — 10i;2 + 26, and since f{m,n) = 
■^1 + "^2 + Vx, we conclude that f{m, n) < AviV2 — 9f 1 — 9v2 + 26. The maximum number of 
facets can be reached when h{£i,£2) vanishes. This occurs when £1 = £2 = 0. That is, when 
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the number of edges of G{P) and G{Q) is maximal. This concludes the proof of the upper 
bound of Theorem 14. II for the special case k = 2. □ 

Corollary 4.4. The maximum number of facets can be attained only when the number of 
edges of each of P and Q is maximal for the given number of facets. 



4.2 The Lower Bound for k = 2 

Given two integers m > 4 and n > 4, we describe how to construct two polytopes in 
with m and n facets respectively, such that the number of facets of their Minkowski sum is 
exactly Amn — 9m — 9n + 26, establishing the lower bound of Theorem 14.11 for the special 
case k = 2. More precisely, given i, we describe how to construct a skeleton of a polytope 
Pi with i facets, 3z — 6 edges, and 2i — 4 vertices, and prove that the number of facets of the 
Minkowski sum of Pm and P„, properly adjusted and oriented, is exactly 4mr;, — 9m — 9r2 + 26. 
As in the previous sections we mainly operate in the dual space of Gaussian maps. However, 
the construction of the desired Gaussian maps described below is an involved task, since not 
every arrangement of arcs of great circles embedded on the unit sphere, the faces of which 
are convex and the edges of which are strictly less than vr long constitutes a valid Gaussian 
map. 

We defer the treatment of the special case z = 4 to the sequel, and 
start with the general case i > 5. The figure to the right depicts the 
Gaussian map of P5. We use the subscript letter i in all notations 
Xi to identify some object X with the polytope Pj. For example, 
we give the Gaussian map G{Pi) of Pj a shorter notation Gj, but in 
this paragraph we omit the subscript letter in all notations for clarity. 
First, we examine the structure of the Gaussian map G of P to better 
understand the structure of P. Let V and E denote the set of vertices 
and edges of G, respectively. Recall that the number of vertices, edges, and faces of G is i, 
3i — 6, and 2i — 4, respectively. The unit sphere, where G is embedded on, is divided by 
the plane y = into two hemispheres C {{x,y, z) \ y < 0} and C {{x,y, z) \ y > 0}. 
Three vertices, namely u, v, and w, lie in the plane x = 0. u is located very close to the 
pole (0,0,-1). It is the only vertex (out of the i vertices) that lies in H~^. v is located 
exactly at the opposite pole (0, 0, 1), and w lies in H~ very close to v. None of the remaining 
i — 3 vertices in \ {u, v, w} lie in the plane x = 0; they are all concentrated near the pole 
(0,0, —1) and lie in H . The edge uv, which is contained in the plane x = 0, is the only 
edge whose interior is entirely contained in H^. Every vertex in \ {u, v, w} is connected 
by two edges to v and w, respectively. These edges together with the edge uw, contained in 
the plane a; = 0, form a set of 2z — 5 edges, denoted as E' = E \ {uv}. The length of each of 
the edges in E' is almost vr, due to the near proximity of u, v, and w to the respective poles. 

It is easy to verify that if the polytope P is not degenerate, namely, its affine hull is 
3-space, then any edge of G{P) is strictly less than vr long. Bearing this in mind, the main 
difficulty in arriving at a tight-bound construction is forcing sufficient edges of the set E' of 
the Gaussian map of one polytope to intersect sufficient edges of the set E' of the Gaussian 
map of the other polytope. The remaining pair of edges, one from each Gaussian map. 
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contributes a single intersection to the total count. As shown below, this is the best one can 
do in terms of intersections. 

The number of facets of the Minkowski sum of Pm and P„ is max- 
imal, when the number of vertices in the overlay of Gm and G„ is 
maximal. This occurs, for example, when one of them is rotated 90° 
about the Y depicted on the left for the case of m = = 5. 

In this configuration, all the 2m — 5 edges in E'^ intersect all the 

^^^^^jfl 2n — 5 edges in E'^. All intersections occur in H~. In addition, the 
edge uVm intersects the edge uVn- The intersection point lies in 
exactly at the pole (0,1,0). Counting all these intersections results 
with (2m — 5)(2n — 5) + 1 = 4mn — 10m — lOn + 26. Adding the 
original vertices of G{Pm) and G{Pn), yields the desired result. 

Next, we explain how Pj, z > 5 is constructed to match the de- 
scription of Gi above. The construction of Pj is guided by a cylinder. 
All the vertices of P, lie on the boundary of a cylinder the axis of 
which coincides with the Z axis. We start with the case i = 5, and show how to generalize 
the construction for i > 5. The special case i = 4 is explained last. 



Figure 4.2: The over- 
lay of G5 and Gg, 
where G5 is G5 rotated 
90° about the Y axis. 



4.2.1 Constructing P5 

Figure shows various views of P5. Recall that P5 has 6 vertices, denoted as fo, f 1, . . . , fs, 
and 9 edges. We omit the subscript digit 5 in all the notations through the rest of this 
subsection for clarity. Let f if 2 • • • f n denote the face defined by the sequence of vertices 
vi,V2, . . . ,Vn on the face boundary. The projection of all vertices onto the plane 2 = lie 
on the unit circle. As a matter of fact, the entire face = VqVivWs lies in the plane z = 0. 
It is mapped under G to the vertex v = G{f""). Similarly, the faces /" = V5V4V2V1 and 
= are mapped under G to the vertices u = G^f^) and w = G{f^), respectively. 

Consider the projection of the vertices onto the plane z = best seen in Figure I4.3r b). 
Once the projection vi^ of ^5 is determined as explained below, vq is placed exactly on the 
bisector of Zv'^ovi. The vertices ^4, V3, and V2 are the reflection of the vertices v^, vo, and 
Vi respectively through the plane x = 0. 

Two parameters govern the exact placement of (and V4). One is the size of the exterior- 
dihedral angle at the edge vqVs, denoted as a, that is, the length of the geodesic-segment 
that is the mapping of the edge vw of G. This angle is best seen in Figure I^Sl fc). Notice, 
that the Z axis is scaled up for clarity, and the angle in practice is much smaller. The other 
parameter is the size of the angle /? = Zv'^ov'^, where v'^ and v'^ are the projections of and 
f5 respectively onto the plane z = 0. This is best seen in Figures [4.3( b) and (e). Given m 
and n, these angles for each of P^ and P„ depend on both m and n. For large values of m 
and n the values of a and (3 should be small. For example, setting a = [3 = 10° is sufficient 
for the case m = n = 5 depicted in Figure 14.21 The actual setting is discussed below after 
the description of the general case i > 5. 




4.2.2 Constructing Pi,i>5 



We construct a polytopc, such that two facets are visible 
when looked at from ^ = 00, and i — 2 facets are visible when 
looked at from z — —00. First, we place the projection of 
all vertices onto the plane z — along the unit circle, and 
denote the projection of a vertex v as v'. The projection of 
the vertices Vj^^, Vj^, Vj^, Vj.^, Vj^, and f^v,, where jo = 0, ji = 
[{t - 2)/2j , j2 = - 2)/2j + 1, j3 = ^ - 2, J4 = L(3^ - 7)/2j , 
and js = [{3i — 7)/2j +1, are placed at the same locations 
as those of the corresponding vertices of P5, as depicted on 

the right. The projection of the remaining vertices are placed on the arcs Vj^,Vjg, Vj^/uJ^^, 




Vj^jVj^, and Vj^, v',^ in cyclic order. 



The angle 7 = ZvjgOVj^-i is another parameter that governs the final configuration of 
Pi. Once the placement of the projection of fji-i is determined, the projections of the 
vertices f jo+i? "^io+Si • • • ! '^ji-2 are arbitrarily spread along the open arc Vj^^TvJ^-i. The vertex 
placement along the arc Vj^,VjQ must be a symmetric reflection of the vertex placement along 
the arc vZ/u7^. This guarantees that all the quadrilateral facets are planar. Similarly, the 




vertex placement along the arc Vj^ , vj^ is a symmetric reflection of the vertex placement along 

the arc Vj.^, Vj^. For large values of m and n the angle 7 should be small as explained below, 
implying that the projection of the vertices are concentrated near Vj^ and Vj,^, (which lie on 
the bisectors of Zvjj^ovj^ and Zvj^ov'j^, respectively). Figure |13] depicts the cases i = 10, and 
i = 11. In these examples we force a regular placement, which is sufficient in many cases. 
As in the case of i = 5, the face = vq, vi, . . . , Vi-2, represented by the vertex Vi of Gi, 
lies in the plane z = 0. The exterior-dihedral angle a at the edge Vj^Vj.^ is made small, so 
that the vertex Wi of Gi representing the adjacent face fl" = Vj^,Vj^+i, . . . , f2j-5, Vq, is kept 
in close proximity to Vi. 

Given m and n, three parameters per polytope listed below govern the final configurations 
of Pm and Pn- 

1. the exterior-dihedral angle a at the edge vj^/v]^, 

2. the angle P = Zvj^oVj^, and 

3. the angle 7 = Zvj^oVj-^^i. 

The settings of these angles must satisfy certain conditions, which in turn enable all the 
necessary intersections of edges in the Gaussian map of the Minkowski sum. We denote the 
face Wjg+ifjgfjjf adjacent to /" by The vertex x = G{f^) is the nearest vertex to u. 
The ^/-coordinate of the vertex w„ must be greater than the y-coordinate of the edge xVm 
at z = in P^'s coordinate system. Similarly, the y-coordinate of the vertex Wm must be 
greater than the y-coordinate of the edge xvn at 2; = in P„'s coordinate system.^ This is 
best seen in Figure 1431 (c). The values of the ^/-coordinates of Wn and Wm are simply sm{an) 
and sin(am), respectively. The value of the ^/-coordinate of the edge xVm at z = however 
depends on all the three parameters am, Pm, and 7^- Similarly, the y-coordinate of the edge 
XV n at 2; = in the respective coordinate system depends on Pn, and 7„. Instead of 
deriving an expression that directly evaluates these ^-coordinates, we suggest an iterative 
procedure that decreases the angles at every iteration until the conditions above are met, 
and argue that this procedure eventually terminates, because at the limit, we are back at 
the case where m = n = 5, for which valid settings exist. 



•^The rotation of, say P„, is performed about the Y axis. Thus, it has no bearing on y-coordinates. 
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Figure 4.5: (a) The Minkowski sum Mn^n = Pu P[^, where is Pu rotated 90° about the Y 
axis, (b) The Gaussian map of Afn^n looked at from z = oo. (c) A scaled up view of the Gaussian 
map of Mil 11 looked at from z = oo. (d) The Gaussian map of Mu^u looked at from y = — oo. 



4.2.3 Constructing P4 




Recall that P4 has 4 facets, 6 edges, and 4 vertices. Therefore, it cannot 
be constructed according to the prescription provided in the previous 
section. Applying the same principles though, we place two vertices 
of G4 near the pole (0, 0,-1), and two vertices near the opposite pole 
(0,0, 1). One edge, which connects a vertex near one pole to a ver- 
tex near the other, lightly shaded in the figure on the left, is entirely 
contained in H~^. The other three edges that connect vertices near 
opposite poles mostly lie in . They form a set of 2z — 5 = 3 edges, denoted as E'^. The 
length of every edge in E'^ is almost vr. In contrast to the case z > 5, two out of the three 
edges in E'^ cross the plane y = 0. Namely, small sections of them lie in H~^. As in the case 
i > 4, one edge, the lightly shaded one, is entirely contained in H~^. 

We construct P4, such that the two facets = V0V1V2 and 

= V0V2V3 are visible when looked at from z = 00, and when 
looked at from z = —00, the remaining two facets = v^viVq 
and = V3V2V1 are visible. As depicted on the right, the 
projection of all four vertices onto the plane 2 = lie on the 
unit circle. The vertices vo and V2 lie on the plane z = 0, and 
the vertices vi and V3 lie in a parallel plane. The distance 
between the planes is small to form small exterior-dihedral 
angles at the edges V0V2 and vTv^- 

As in the general case, two parameters govern the exact placement of Vi, V2, and ^3. One 
is the size of the exterior-dihedral angle at the edge VoV2- The other parameter is the size of 
the angle Zv20v'^, where ^3 is the projection of V3 onto the plane z = 0. The sizes of these 
angles are determined by the same rationale as in the general case. 

This concludes the proof of the lower bound of Theorem 14.11 for the special case k = 2. 
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4.3 Maximum Complexity of Minkowski Sums of Many 
Polytopes 

Let Pi, P2, . . . , Pk be a set of k polytopes in R^, such that the number of facets of Pi is rrii 
for i = 1,2, . . . , k. In this section we present a tight bound on the number of facets of the 
Minkowski sum M = Pi (B P2 (B ■ ■ ■ (B Pk generahzing the arguments presented above for 
k = 2. 



4.3.1 The Lower Bound 




such that 



> 4, 



we 



de- 



Given k positive integers mi,m2, ■ ■ ■ ,mk, 
scribe how to construct k polytopes in with corresponding number 
of facets, such that the number of facets of their Minkowski sum is 
,(2mj — 5)(2mj — 5) + (2) + Yli=i '^i- More precisely, 



exactly 



<i<j<k^ 



given i, we describe how to construct a skeleton of a polytope Pj with 
i facets, 3i — 6 edges, and 2z — 4 vertices, and prove that the number of 
facets of the Minkowski sum M = Pi® P2® ■ ■ ■ ® Pk oi the k polytopes 
properly adjusted and oriented is exactly the expression above. We use 
the same construction described in Section I4.2[ 

The number of facets in the Minkowski sum of Pi,i = 1,2, ... ,k 
is maximal, when the number of vertices in the overlay of Gi,i = 
1,2, ... ,k is maximal. This occurs, for example, when Gi is rotated 
180° (i — l)/k about the Y axis for z = 1, 2, . . . , fc, as depicted on the 
left for the case of mi = 1712 = = 4. (Recall, that E'^ refers to 
the set of edges that span the lowest hemispheres, and its cardinality 
is smaller than cardinality of E by one.) In this configuration, all the 
5 edges in E^ intersect all the 2mj — 5 edges in Ej, for 1 < i < j < k. All intersections 
occur in H~. In addition, the edge uvrm intersects the edge uVmj ioi 1 < i < j < k. These 
intersection points lie in near the pole (0, 1,0). Counting all these intersections results 
with X]i<i<i<fc(2"^j "~ 5)(2mj — 5) + (2) . Adding the original vertices of G{Pi), 
yields the bound asserted in Theorem 14. 1[ 



Figure 4.6: The 
overlay of the Gaus- 
sian maps of three 
tetrahedra rotated 
about the Y axis 
0°, 60°, and 120°, 
respectively. 



2m,- 



1,2, 



, k. 



4.3.2 The Upper Bound 



We can apply the special case A; = 2 of Theorem 14.11 to obtain 

/(mi, ma, ... , m^) < /(mi, /(mg, mg, . . . , m^)) 

< 4mi/(m2, ma, ... , ruk) - 9mi - 9/(m2, mg, 

k 

<4'= J]m, + ... 

i=l 



.mk) + 26 



However, we can apply a technique similar to the one used in Section 14.11 and improve this 
upper bound, but first we must extend Lemma IT3l 
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Lemma 4.5. Let Gi, G2, . . . ,Gk be a set of k Gaussian maps of convex polytopes, and let 
G be their overlay. Let fi denote the number of faces of Gi, and let f denote the number of 
faces of G. Then, f < Ei<.<,<. f^f^ ' (k - 2) Ei<,<, /. + (fc - - 2). 

Proof. Let us choose two antipodal points on the sphere S^, such that no arc of the overlay 
is aligned with them. In particular, the points are in the interior of two distinct faces. We 
consider these two points to be the north pole and south pole of the sphere, and define the 
direction west to be the clockwise direction when looking from the north pole toward the 
south pole. We define a western-most corner to be a pair of a face and one of its vertices, 
which is to the west of all of its other vertices. Apart from the two faces, which contain the 
poles, any face has a unique western-most corner, since no edge is aligned with the poles, 
and all faces of any Gaussian map are spherically convex. So for any overlay with / faces, 
there are / — 2 such western-most corners. 

The maximal number of faces is attained when the overlay G is non-degenerate. Thus, a 
vertex of G is either the intersection of two edges of some distinct Gi and Gj, or a vertex of 
some Gj. Therefore, a western- most corner for a face of G is either a western- most corner for 
the overlay of some Gi and Gj, or a western-most corner for some Gi, in which case it also 
is a western-most corner for any overlay involving Gi. The number of western- most corners 
in the Gaussian map Gi is fi — 2, and the maximal number of western-most corners in the 
overlay of some Gi and Gj is fifj — 2. 

We can therefore write: 

k 

f< Yl (M-2)-(fc-2) J](/,-2) + 2. 

l<i<j<k 1=1 

This corresponds to summing the western-most corners appearing in the overlay of all pairs 
of Gaussian maps, and subtracting [k — 2) times the western-most corners appearing in all 
original Gaussian maps, since each of them appeared [k — 1) times in the first sum. Finally, 
we have: 

k k 

J2 (M-2)-(fc-2) J](/,-2) + 2= hf,~{k-2)Y,h + {k-l){k-2) . 

l<i<j<k 1=1 l<i<j<k i=l 

□ 

Let Pi, P2, ■ ■ ■ ,Pk be k polytopes in M'^ with mi, m2, . . . , facets, respectively. Let 
G{Pi) denote the Gaussian map of Pi, and let Vi, Ci, and fi denote the number of vertices, 
edges, and faces of G{Pi), respectively. Let denote the number of intersections of edges 
of G{Pi) and edges of G{Pj), i ^ j m. G{M). Applying the same technique as in Section WA] 
that is, counting the total degrees of vertices in G{M) implies that ^j=i Cj + 2t>^ = e. Using 
Euler's formula, we get Yl^=i + = f + v — 2. Applying Lemma 14.51 and respecting 
V = X^KKfc"^* + "^a;' '^^^ bound from above: 



k k 

v.< Yl f^fJ-ik-2)Yf^ + ik-l){k-2) + Yi^,-ei)-2 . (4.2) 

l<i<j<k 2=1 i=l 



76 



Chapter 4- Exact Complexity of Minkowski Sums 



According to Corollary 14.41 the maximum number of facets of the Minkowski sum of two 
polytopes is attained when the number of edges of each summand is maximal. We need to 
establish a similar property for the general case. Generalizing the derivation procedure in 
Section H?T| we introduce k non-negative integers = 1,2, k, such that Cj = 3f « — 6 — ii 
and fi = 2vi — 4: — ii. Substituting Cj in (14.21) we get: 

k k 

v,< Yl I\fj-ik-2)Y,J\ + ik-l)ik-2) + Y,ivi-^Vi + 6 + ii)-2 

l<i<j<k 1=1 i=l 

k k k 

< E M-(^-2)$^/.-$^(2i;.-5) + 5^£, + fc2-2fc. (4.3) 

l<i<jr<A; i=l i=l i=l 

Substituting /j in (14.31) we get: 

k k 

Vx< {2vi-4:-^i){2vj-A~^j)-{k-2)Y{2vi-A-^,)-Y{2v,~b-^i)+k'^ -2k 

l<i<j<k 1=1 i=l 



J2 {2v.-5){2v,-5)+Q+h{i^ 



, ' • ' , ' 



l<i<j<k 

where 

h{£i, 4, . . . , 4) = Yl - ^^-(^^^ - 5) - - 5)) 

k 

= E^^(E(V2-(2t;,-5))). 

i=l j^i 

Connectivity of G{Pi) implies that £i < 2vi — 5, which in turn implies that h{ii,i2, • • • , 4) < 
0. Thus, we have: 



Yi i'2v,-5)i2vJ-5)+^^^ . 

l<i<j<k ^ ^ 



(4.4) 



We conclude that the exact maximum number of facets of the Minkowski sum of k 
polytopes cannot exceed ^i<j<j<fc(2mi — 5)(2mj — 5) + ^^<j<^mj+ (2), which completes the 
proof of Theorem 14. 1[ For example, the exact maximum number of facets of the Minkowski 
sum of k tetrahedra is 5k'^ — k. For k = 2 the expression evaluates to 18. 



Divide each difEculty into 
as many parts as is feasible 
and necessary to resolve it. 



Rene Descartes 
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Assembly Planning 



Assembly partitioning with an infinite translation is the application of an infinite translation 
to partition an assembled product into two complementing subsets of parts, referred to as 
subassemblies, each treated as a rigid body. We present an exact implementation of an 
efficient algorithm based on the general framework described in jHLWOO| IWL94] to obtain 
such a motion and subassemblies given an assembly of polyhedra in M^. As usual, we do not 
assume general position. Namely, we handle degenerate input, and produce exact results. As 
often occurs, motions that partition a given assembly or subassembly might be isolated in the 
infinite space of motions. Any perturbation of the input or of intermediate results, caused 
by, for example, imprecision, might result with dismissal of valid partitioning-motions. In 
the extreme case, there is only a finite number of valid partitioning motions, as occurs in the 
assembly shown in Figure 15.11 no motion may be found, even though such exists. Proper 
handling requires significant enhancements applied to the original algorithmic framework. 
The implementation is based on software components introduced in Chapters [2] and [31 They 



Figure 5.1: (a) The Split Star assembly, and (b),(c), and (d) the Split Star six parts divided into three 
pairs of symmetric parts. The six parts are named according to their color R{ed), G(reen), B{\ue), 
P(urple), y(ellow), and r(urquoise). 




(a) 





(d) 
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paved the way to a complete, efficient, and concise implementation. 



5.1 Introduction 

Assembly planning is the problem of finding a sequence of motions that transform the ini- 
tially separated parts of an assembly to form the assembled product. The reverse order of 
sequenced motions separates an assembled product to its parts. Thus, for rigid parts, assem- 
bly planning and disassembly planning refer to the same problem, and used interchangeably. 
In this chapter we concentrate on the case where the assembly consists of polyhedra in 
and the motions are infinite translations. 

The motion space is the space of possible motions that assembly parts may undergo. For 
each motion in a motion space, a subset of parts of a given assembly may collide with a 
different subassembly, when transformed as a rigid body according to the motion. Pairs of 
subassemblies that collide constitute constraints. The motion space approach dictates the 
precomputation of a decomposition of a motion space into regions, such that the constraints 
among the parts in the assembly are the same for all the motions in the same region. All 
constraints over a region are represented by a graph, called the directional blocking graph 
(DBG) [ WL94j . The collection of all regions in a motion space together with their associated 
DBGs, collectively referred to as the nondirectional blocking graph (NDBG), can be used to 
obtain assembly (or disassembly) sequences. 

The general framework and some of the techniques presented here have already been 
described in a series of papers and reports published in the past mainly during the late 90 's. 
Halperin, Latombe, and Wilson made the connection between previously presented tech- 
niques that had used the motion space approach, and introduced a unified general frame- 
work [HLWnnj at the end of the previous millennium. Only few publications related to this 
topic appeared ever since, to the best of our knowledge, which creates a long gap in the time 
line of the respective research. We certainly hope that the tools exposed in this chapter will 
help revive the research on algorithmic assembly planing, a research subject of considerable 
importance. Moreover, we believe that the machinery presented here, together with other 
recent advances in the practice of computational-geometry algorithms, can more generally 
support the development of new and improved techniques in algorithmic automation [2J. 

Solution to the assembly planning problem enables better feedback to designers. It pro- 
vides a design team with additional tools to asses a design, prior to the construction of 
physical mock-ups, and helps creating products that are more cost-effective to manufacture 
and maintain. This is emphasized in light of the strategy to "plan anywhere, build any- 
where" many Computer Aided Design (CAD) and Computer Aided Manufacturing (CAM) 
companies are trying to adopt. Assembly sequences are also useful for selecting assembly 
tools and equipment, and for laying out manufacturing facilities. 

We restrict ourselves to two-handed partitioning steps, meaning that we partition the 
given assembly into two complementing subsets each treated as a rigid body. Even for two- 
handed partitioning, if we allow arbitrary translational motions (and not restrict ourselves 
to infinite translations) the problem is NP-hard [ KK95] . The more general problem of 
assembly sequencing, namely planning a total ordering of assembly operations that merge 
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the individual parts into the assembled product, is PSPACE-hard, even when each part is a 
polygon with a constant number of vertices [Nat88j . 

Notice that the problem that we address in this chapter, namely partitioning with infinite 
translations, is technically considerably more complex than partitioning with infinitesimal 
motions. Although the latter may sound more general, as it handles infinitesimal transla- 
tions and rotations, it is far simpler to implement, since it deals only with constraints that 
can be described linearly. Thus, the problem can be reduced to solving linear programs. 
Indeed, there are several implementations for partitioning with infinitesimal motions (see, 
e.g., iGHH+QSi ISS94j ). but none that we are aware of dealing robustly with infinite trans- 
lations. The shortcoming of using infinitesimal motions only is that suggested disassembly 
moves may not be extendible to practical finite-length separation motions. 

Infinite-translation partitioning was not fully robustly implemented until recently, in spite 
of being more useful than infinitesimal partitioning, most probably due to the hardship of 
accurately constructing the underlying geometric primitives. What enables the solution that 
we present here, is the significant headway in the development of computational-geometry 
software over the past decade, the availability of stable code in the form of the Computa- 
tional Geometry Algorithms Library (Cgal) in general and code for Minkowski sums and 
arrangements in particular [WFZHOTbj . 

The implementation presented in this chapter is based on the Arrangement_on_surf ace_2 
package of Cgal; see Chapter [2] for more details. The implementation uses in particular 
arrangements of geodesic arcs embedded on the sphere; see Section 12.61 The ability to ro- 
bustly construct such arrangements, and carry out exact operations on them using only 
(exact) rational arithmetic is a key property that enables an efficient implementation. The 
implementation exploits supported operations, and requires additional operations, e.g., cen- 
tral projection of polyhedra, which we implemented. We plan to make these new components 
available as part of a future public release of Cgal as well. 



5.1.1 Split Star Puzzle 

We use the assembly depicted in Figure 15.11 as a running example throughout the chapter. 
The name "Split Star" was given to this shape by Stewart Coffin in one of his Puzzle Craft 
booklet editions back in 1985. He uses the term puzzle to refer to any sort of geometric 
recreation having pieces that come apart and fit back together. We use it as an assembly. 
He describes how to produce the actual pieces out of wood |Cof06] , and suggests that they are 
made very accurately. He observes that finding the solution requires dexterity and patience, 
when the pieces are accurately made with a tight fit. Even though the assembly seems 
relatively simple, this should come as little surprise, since the first partitioning motion is 
one out of only eight possible translations of four symmetric pairs of motions in opposite 
directions associated with two complementing subassemblies of three parts each. Evidently, 
any automatic process that introduces even the slightest error along the way, will most likely 
fail to compute the correct first motion in the sequence, and falsely claim that the assembly 
is interlocked. 



80 



Chapter 5. Assembly Planning 





The Split Star assembly has the assembled shape of the 
first stellation of the rhombic dodecahedron |Luk57] . illus- 
trated atop the right pedestal in M. C. Escher's Waterfall 
woodcut |Boo82j . Its orthogonal projection along one of 
its fourfold axes of symmetry is a square, while the Star 
of David is obtained when it is projected along one of its 
threefold axes of symmetry, as seen on the left; for more de- 
tails see |Cof06j . The assembly is a space- filling solid when 
assembled. It consists of six identical concave parts. Each 
part can be decomposed into eight identical tetrahedra yielding 48 tetrahedra in total. As 
manufacturing the pieces requires extreme precision, it is suggested to produce the 48 iden- 
tical pieces and glue them as necessary. Each part can also be decomposed into three convex 
polytopes — two square pyramids and one octahedron, yielding 18 polytopes in total. The 
partitioning described in this chapter requires the decomposition of the parts into convex 
pieces. The choice of decomposition may have a drastic impact on the time consumption of 
the entire process, as observed in a different study in |AFH02] , and shown by experiments 
in Section [531 



5.1.2 Chapter Outline 

The rest of this chapter is organized as follows. The partitioning algorithm is described 
in Section 15.21 along with the necessary terms and definitions. In Section 15.31 we provide 
implementation details. Section 15.41 presents optimizations that are not discussed in the 
preceding sections, some of which we have already implemented and proved to be useful. We 
report on experimental results in Section 15.51 



5.2 The Partitioning Algorithm 

The main problem we address in this chapter, namely, polyhedral assembly partitioning with 
infinite translations, is formally defined as follows: Let A = {Pi, P2, . . . , Pn} be a set of n 
pairwise interior disjoint polyhedra in M^. A denotes the assembly that we aim to partition. 
We look for a proper subset S G A and a direction d in M^, such that S can be moved as 
a rigid body to infinity along d without colliding with the rest of the parts of the assembly 
A \ 5*. (We allow sliding motion of one part over the other. We disallow the intersection of 
the interior of two polyhedra.) 

We follow the NDBG approach |WL94j . and describe it here using the general formulation 
and notation of [HLWOO J . The motion space in our case, namely the space of all possible 
partitioning directions, is represented by the unit sphere S^. Every point p on defines the 
direction from the center of §^ towards p. Every direction d is associated with the directed 
graph DBG{d) = (V, E) that encodes the blocking relations between the parts in A when 
moved along d as follows: The nodes in V correspond to polyhedra in A; we denote a node 
corresponding to the polyhedron Pi by v{Pi) G V. There is an edge directed from v{Pi) 
to v{Pj), denoted e{Pi,Pj) e E, if and only if the interior of the polyhedron Pj intersects 
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the interior of the polyhedron Pj when Pi is moved to infinity along the direction d, and Pj 
remains stationary. 

The key idea behind the NDBG approach is that in problems such as ours, where the 
number of parts is finite, and any allowable partitioning motion can be described by a small 
number of parameters, there is only a relatively small (polynomial) number of distinct DBGs 
that need to be constructed in order to detect a possible partitioning direction. Stated 
differently, the motion space can be represented by an arrangement comprising a finite 
number of cells each assigned with a fixed DBG. Once this arrangement is constructed, we 
construct the DBG over each cell of the arrangement, and check it for strong connectivity. 
A DBG that is not strongly connected is associated with a direction, or a set of directions in 
case the cell is not a singular point, that partition the given assembly. The desired movable 
subset S G A is a byproduct of the algorithm that checks for strong connectivity. If all the 
DBGs over all the cells of the arrangement are strongly connected, we conclude that the 
assembly is interlocked, as a subset of the parts in A that can be separated from the rest of 
the assembly by an infinite translation does not exist. 

Next we show how to construct the motion-space arrangement and compute the DBG 
over each one of the arrangement cells. Each ordered pair of distinct polyhedra < Pi, Pj > 
defines a region Qij on which is the union of all the directions d, such that when Pj is 
moved along d its interior will intersect the interior of Pj. How can we effectively compute 
this region? Let Mij denote the Minkowski sum Pj Q) {—Pi) = {b — a\a G Pi, b G Pj}. We 
claim that the central projection of Mij onto is exactly Qij. 

Lemma 5.1. A direction d is in the interior of the central projection of Mij onto §^ if and 
only if when Pi is moved along d its interior will intersect the interior of Pj. 

Proof. Let d be some direction in the central projection of Mij onto S^. In other words, there 
exists a point m G Mij, such that m = s ■ d, for some positive scalar s. As m is in Mij, there 
exist two points Pi G Pi and Pj G Pj, such that m = pj — pi. Thus, Pj = Pi + s ■ d, meaning 
that the point Pi intersects pj when moved along d. A similar argument can be used to show 
the converse. □ 

Next, we describe how, given two polyhedra Pi and Pj, we compute the region Qij, using 
robust and efficient hierarchy of building blocks, which we have developed in recent years. 
The existing tools that we use are (i) computing the arrangement of spherical polygons 
[BFH+OTl [FSHOSbl IFSH08al IWFZH07bj . and (ii) construction of Minkowski sums of convex 
polytopes |BFH+n7i IFHOTl IFSHOSbj IFSHOSaj . We also need some extra machinery, as 
explained below. 

Assume Pi is given as the union of a collection of (not necessarily disjoint) convex 
polytopes PI, P2, . . . , P^., and similarly Pj is given as the collection of convex polytopes 
Pi, Pi, ...,Pi. It is easily verified that M,, = {jk=i,...,ma=i,-,m, Pi ® (-^fc)- So we com- 
pute the Minkowski sum of each pair Pi © {—PI), and centrally project it onto S^. Finally, 
we take the union of all these projections to yield Qij. 

There are several ways to effectively compute the central projection of a convex polyhe- 
dron C (one of the polytopes M^-^ = Pj © {—PI)) from the origin onto S^. We opted for the 
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following. An edge e of C is called a silhouette edge, if the plane vr through the origin and e 
is tangent to C at e; namely, vr intersects C in e only. We assume for now that no tangent 
plane contains a facet of C; we relax this assumption in Section [5.3.5^ where we provide a 
detailed description of the procedure. We traverse the edges of C till we find a silhouette 
edge ei. One can verify that the silhouette edges form a cycle on C. We start with ei, and 
search for a silhouette edge adjacent to ei. We proceed in the same manner, till we end up 
discovering ei again. Projecting this cycle of edges onto §^ is straightforward. 

All the boundaries of the regions Qij form an arrangement of geodesic arcs on the sphere. 
We traverse the motion-space arrangement in say a breadth-first fashion. For the first face 
we check which ones of the regions Qij contain it. We construct the corresponding DBG 
and check it for strong connectivity. If it is not strongly connected, we stop and announce 
a solution as described above. Otherwise we move to an adjacent feature of the current 
face. During this move we may have stepped out from a set of regions Qij, and may have 
stepped into a new set of regions Qij. We update the current DBG according to the regions 
we left or entered, test the new DBG for strong connectivity, and so on till the traversal of 
all the arrangement cells is completed. Notice that it is important to visit also vertices and 
edges of the arrangement, since the solution may not lie in the interior of a face. Indeed, in 
our Split Star example, solutions are on vertices of the arrangement. Without careful exact 
constructions, such solutions could easily be missed. 

5.3 Implementation Details 

The implementation of the assembly-partitioning operation consists of eight phases listed 
below. They all exploit arrangements of geodesic arcs embedded on the sphere |BFH+07 
IFSH08aj in various ways. The Arrangement_on_surf ace_2 package of Cgal already sup- 
ports the construction and maintenance of such arrangements, the computation of union of 
faces of such arrangements, the construction of Gaussian maps of polyhedra represented by 
such arrangements, and the computation of their Minkowski sums. It provides the ground 
for efficient and generic implementation of the remaining required operations, such as central 
projection. 

1. Convex Decomposition 

2. Sub-part Gaussian map construction 

3. Sub-part Gaussian map reflection 

4. Pairwise sub-part Minkowski sum construction 

5. Pairwise sub-part Minkowski sum projection 

6. Pairwise Minkowski sum projection 

7. Motion-space construction 

8. Motion-space processing 
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Figure 5.2: The Split Star six parts decomposed into three convex sub-parts each. 

The partitioning process is implemented as a free function that accepts as input an 
ordered hst of polyhedra in M^, which are the parts of the assembly. Each part is represented 
as a polyhedral mesh in M^; see Section [3.2.11 for a definition. We proceed with a detailed 
discussion of the implementation of each phase. 

We deal below with various details that are typically ignored in reports on geometric 
algorithms (for example, under the general position assumption). However, in assembly 
planning, or more generally in movable-separability problems |Tou85j in tight scenarios, 
much of the difficulty shifts exactly to these technical details and in particular to handling 
degeneracies. This is especially emphasized in Phases 5 and 6 (Subsections 15.3.51 and 15.3.61 
respectively), but prevails throughout the entire section. 

5.3.1 Convex Decomposition 

We decompose each concave part into convex polyhedra referred to as sub-parts. The output 
of this phase is an ordered list of parts, where each part is an ordered list of convex sub- 
parts represented as polyhedral surfaces. Each polyhedral surface is maintained as a Cgal 
Polyhedron_3 |Ket07aj data structure, which consists of vertices, edges, and facets and 
incidence relations on them |Ket99] . A part that is convex to start with is simply converted 
into an object of type Polyhedron_3. 

A new package of Cgal that supports convex decomposition of polyhedra has been 
recently introduced [Hac07j . but has not become publicly available yet. As we aim for a 
fully automatic process, we intend to exploit such components, once they become available, 
and study their impact. For the time being we resorted to a manual procedure. A simple 
decomposition of the Split Star parts used in the running example is illustrated in Figure 15^ 

5.3.2 Sub-part Gaussian Map Construction 

We convert each sub-part represented as a polyhedral surface into a Gaussian map repre- 
sented as an arrangement of geodesic arcs embedded on the sphere, where each face / of the 
arrangement is extended with the coordinates of its associated primal vertex v = G~^{f), 
resulting with a unique representation; see Section 13.11 for the exact definition of Gaussian 
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Figure 5.3: Samples of the Gaussian maps of sub-parts of the Split Star assembly. The bottom row 
contains the reflections of the Gaussian maps at the top row. 

maps and see Section 13.2! the exact procedure to construct one from a polytope. 

The output of this phase is an ordered hst of parts, where each part is an ordered hst of 
the Gaussian maps of the convex sub-parts. Figure 15.3! depicts the Gaussian maps of six of 
the 18 polytopes that comprise the set of sub-parts of our Spht Star assembly. 



5.3.3 Sub-part Gaussian Map Reflection 

We reflect each sub-part through the origin to obtain — P^. This operation can be per- 
formed directly on the output of the previous phase, reflecting the underlying arrangements 
of geodesic arcs embedded on the sphere, which represent the Gaussian maps, while handling 
the additional data attached to the arrangement faces. As a matter of fact, this phase can 
be reduced as part of an optimization discussed in Section 15.41 

The output of this phase is an ordered list of parts, where each part is an ordered list of 
Gaussian maps of the reflected convex sub-parts. Figure 15.3! depicts the Gaussian maps of 
six of the 18 polytopes that comprise the set of reflected sub-parts of the Split Star example. 



5.3.4 Pairwise Sub-part Minkowski Sum Construction 



We compute the Minkowski sums of the pairwise 
sub-parts and reflected sub-parts. Aiming for an ef- 
ficient output sensitive algorithm, the construction 
of an individual Minkowski sum from two Gaussian 
maps represented as two arrangements respectively 
is performed by overlaying the two arrangements. 
When the overlay operation progresses, new vertices, 
edges, and faces of the resulting arrangement are cre- 
ated based on features of the two operands. When a 
new feature is created its attributes are updated. There are ten cases that must be handled; 
see Sections 13. 2. 2l for details. The Arrangement_on_surf ace_2 package conveniently supports 
the overlay operation allowing users to provide their own version of these ten operations. The 
overlay operation is exploited below in several different variants of arrangements of geodesic 



Construct Pairwise Sub-part 

Minkowski Sums 
for i in {1, 2, . . . , n} 
for j in {1,2,..., n} 
if i == j continue 
for k in {1, 2, . . . , rrii} 

for £ in {1,2,..., rrij} 
'-PL) 
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Figure 5.4: Samples of the pairwise Minkowski sums of sub-parts of the Split Star assembly. The 
middle row contains six Minkowski sums. The top row contains the corresponding Gaussian maps. The 
bottom row contains the corresponding central projection of the Minkowski sums on S^. 

arcs embedded on the sphere. Each apphcation requires the provision of a different set of 
those ten operations. 

The output of this phase is a map from ordered pairs of distinct indices into hsts of 
Minkowski sums represented as Gaussian maps. Each ordered pair <i,j>,i ^ j is associated 
with the hst of Minkowski sums {M^ \ k = 1,2,..., rrii, i = 1,2, ... ,mj}. In case of our Spht 
Star the map consists of 30 entries that correspond to all configurations of ordered distinct 
pairs of parts. Each entry consists of a list of nine Minkowski sums, that is, 270 Minkowski 
sums in total. 



5.3.5 Pairwise Sub-part Minkowski Sum Projection 



Project Pairwise Sub-part ~ centrally project all pairwise sub-part Minkowski 

Minkowski sums ^^^^ computed in the previous phase onto the 
sphere. Each projection is represented as an arrange- 
ment of geodesic arcs on the sphere, where each cell 
c of the arrangement is extended with a Boolean flag 
that indicates whether all infinite rays emanating 
from the origin in all directions c? G c pierce the cor- 
responding Minkowski sum. As the Minkowski sums 
are convex, their spherical projections are spherically 



for z in {1, 2, ... , n} 
for j in {1, 2, ... , n} 
if i == j continue 
for k in {1, 2, . . . , rrii} 

for £ in {1, 2, ... , mj} 



QZ = WO}ect{Ml^,) 



convex. 



Given a convex Minkowski sum C, we distinguish between four different cases as follows: 

1. The origin is contained in the interior of a facet of C. 

2. The origin lies in the interior of an edge of C. 

3. The origin coincides with a vertex of C. 
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4. The origin is separated from C. 

Computing the projection of a convex polytope C can be done efficiently using dedicated 
procedures that handle the four cases, respectively. Recall that C is represented as a Gaussian 
map, which is internally represented as an arrangement of geodesic arcs embedded on the 
sphere. We traverse the vertices of the arrangement. For each vertex v we extract its 
associated primal facet / = G~^{y). We dispatch the appropriate computation based on the 
relative position of the origin with respect to the supporting plane to /, and the supporting 
plane to adjacent facets of /. 

If the origin is contained in the interior of a facet / of C, the projection of 
the silhouette of C is a great (full) circle that divides the sphere into two hemispheres. The 
normal to the plane that contains the great circle is identical to the normal to the supporting 
plane to /, easily extracted from the arrangement representing the Gaussian map of C . The 
Arrangement_on_surf ace_2 package conveniently supports the insertion of a great circle, 
provided by the normal to the plane that contains it, into an arrangement of geodesic arcs 
embedded on the sphere. 

We omit the implementation details of the following two cases, and proceed to the general 
case. If the origin is separated from C, we traverse all edges of C until we find a silhouette 
edge characterized as follows: Let Vg and Vt be the source and target vertices of some edge e in 
the arrangement representing the Gaussian map of C, and let fg = G~^{vs) and ft = G~^{vt) 
be their associated primal facets, respectively, e is a silhouette edge, if and only if, the origin 
is not in the negative side of the supporting plane to fs and not in the positive side of the 
supporting plane to ft- We start with the first silhouette edge we find, and search for an 
adjacent silhouette edge in a loop, until we rediscover the first one. We 
project only the target vertices of significant silhouette edges, and con- 
nect consecutive projections using arcs of great circle. Let e and e' be 
adjacent silhouette edges. We skip e, if the projections of e and e' lie on 
the same great circle. For example, all but the last adjacent silhouette 
edges incident to a facet supported by a plane that contains the origin are 
redundant, as illustrated in the figure above. Here we skip cq, ei, and 62, 
and project the target vertex of 63. 

The output of this phase is a map from ordered pairs of distinct indices into lists of 
arrangements as described above. Each ordered pair <i,j>,i 7^ j is associated with the list 
of central projections of the pairwise Minkowski sums of Pj's sub-parts and the refiection 
through the origin of Pj's sub-parts. 

5.3.6 Pairwise Minkowski Sum Projection 

For each pair of distinct parts Pj and Pj we compute the union of projections of the pairwise 
Minkowski sums of all sub-parts of part Pj and refiections of all sub-parts of part Pj. 

The output of this phase is a map from ordered pairs of distinct indices into arrangements. 
Each ordered pair < z, j >, i 7^ j is associated with a single arrangement extended as described 
above, that represents the central projection Qij of Mij = Pj © (— Pj)- 
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Minkowski sums Projections 
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for i in {1, 2, ... , n} 
for j in {1,2, . . . ,n} 
if i == j continue 

Qij = 

for k in {1, 2, . . . , rrii} 
for £ in {1,2,..., rrij} 

Qij 



Qij U Qj^^ 



lay of a face in some projection Ql 



We exploit the overlay operation in this phase the 
second time throughout this process, this time in a 
loop. Given two distinct parts Pi and Pj we traverse 
all projections in the set {Q^j^g \ k = 1,2, ... ,mi,i = 
1,2, ... ,mj}, and accumulate the result in the ar- 
rangement Qij . As mentioned in Section I5.3.4[ when 
the overlay operation progresses, new vertices, edges, 
and faces of the resulting arrangement are created. 
When a new face / is created as a result of the over- 
and a face in the accumulating arrangement, the 
Boolean flag associated with /, which indicates whether all directions d E f pierce Mij, is 
turned on, if d pierces M^-^, that is, if the flag associated with the face of g is on. 

The intermediate result of this step are arrangements with potentially redundant edges 
and vertices. It is desired (but not necessary) to remove these cells, as it reduces the time 
consumption of the succeeding operations, which is directly related to the complexity of the 
arrangements. It has even a larger impact when the optimization described in Section [53] is 
applied, as the optimization decreases the number of preceding operation at the account of 
slightly increasing the number of succeeding operations. We remove all edges and vertices 
that are in the interior of the projection, that is all marked edges and vertices. We also 
remove spherically collinear vertices on the boundary of the projection, the degree of which 
decreased below three, as a result of the redundant-edge removal. 






(d) 



(e) 



(f) 



Figure 5.5: Peg-in-the-hole Minkowski sum projections, (a), (b), (c), (d), and (e) are the sub-part 
projection, (f) is the union of the former. 



Cgal also supports Boolean operations applied to general polygons^ and in particular 
the union operation. However, it consumes and produces regularized general polygons; see 
Section 12.7.11 This regularization operation is harmful in the realm of assembly planning. 
Therefore, we work directly on the cells of the arrangements Qij. 
Quite often the projection contains isolated vertices and edges, 
as occurs in the example depicted on the right, referred to as 
"peg-in-the-hole". Here the assembled product is translucently 
viewed from two opposite directions. The blue part is stationary 
and is decomposed into five sub-parts. Figure 15.51 illustrates the corresponding five pairwise 
Minkowski sum projections, and their union. The complement of the union consists of a 
single isolated vertex. 



^The generic code supports point sets bounded by algebraic curves embedded on parametric surfaces 
referred to as general polygons. 
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Recalling our Split Star assembly, the projection of the Minkowski 
sum of the red part and the reflection of the blue part, and its 
reflection, that is, the projection of the Minkowski sum of the blue 
part and the reflection of the red part are depicted on the left. 



B 



'-R) R®{-B) 



5.3.7 Motion-Space Construction 

We compute a single arrangement that represents the motion space, 
where each cell c of the arrangement is extended with a DBG. We 
use the adjacency-matrix storage format provided by Boost [3] to 
represent each DBG. Recall that for a graph with n vertices such 
as ours, an n x n matrix is used, where each element a?- of a DBG associated with cell c is 
a Boolean flag that indicates whether part Pj collides with part Pj when moved along any 
direction d & c. In particular we use the adjacency_matrix class. It implements the Boost 
Graph Library (BGL) |SLL02j interface, which supports, among the other, easy insertions 
of new edges into existing graphs. Handling large assemblies with sparse blocking relations 
may require different representations of DBGs to reduce memory consumption. 

We exploit the overlay operation in this phase the third time similar to its application 
in Section 15.3.61 We traverse all central projections in the set {Qij \l<i<j<n}, 
and accumulate the result in the final motion-space arrangement. As mentioned above in 
Section I5.3.4[ when the overlay operation progresses, new vertices, edges, and faces of the 
resulting arrangement are created. When a new cell c is created as a 
result of the overlay of a face g in some projection Qij, and a cell in 
the accumulating arrangement, the DBG associated with c is updated. 
That is, if the flag associated with g is turned on, we insert an edge 
between vertex i and vertex j into the DBG associated with c. 

Depicted on the left is the motion-space arrangement computed by 
our program for the Split Star assembly. 




5.3.8 Mot ion- Space Processing 



Table 5.1: The Split Star valid 
partitioning directions and cor- 
responding subassemblies. 





Direction 


Subset 


1. 


-1,-1,-1 


rBT 


2. 


-1,-1, 1 


RBT 


3. 


-1, 1,-1 


CrPT 


4. 


-1, 1, 1 


RPT 


5. 


1,-1,-1 


GBY 


6. 


1,-1, 1 


RBY 


7. 


1, 1,-1 


GPY 


8. 


1, 1, 1 


RPY 



We traverse all vertices, edges, and faces of the motion- 
space arrangement in this order, and test the DBG associated 
with each cell for strong connectivity using the Boost global 
function strong_components(). This function computes the 
strongly connected components of a directed graph using Tar- 
jan's algorithm based on depth-first search (DFS) (Tar 72] . The 
set of constraints associated with a vertex f is a proper sub- 
set of the constraints associated with the edges incident to v. 
Similarly, the set of constraints associated with an edge e is a 
proper subset of the constraints associated with the two faces 
incident to e. Therefore, if the DBGs of all vertices are strongly 
connected, we terminate with the conclusion that the 
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assembly is interlocked. Similarly, if we are interested in finding all 
solutions, and the DBGs of all edges are strongly connected, we termi- 
nate, as no further solutions on faces exist. 

For the Split Star assembly, our program successfully identifies all 
the eight partitioning directions depicted on the right along with the 
corresponding subset of parts listed in Table 15.11 

5.4 Additional Optimization 

The reflection of the sub-parts through the origin as described in Section 15.3.31 has been 
naively implemented. The computation is applied to the polyhedral-mesh representation of 
each sub-part. An immediate optimization calls for an application of the reflection operation 
directly on the arrangements that represent the Gaussian map. We are planning to intro- 
duce a generic implementation of the reflection operation that operates on any applicable 
arrangement. This operation alters the incidence relations between the arrangement features 
and their geometric embeddings. For each vertex, it negates its associated point, and inverts 
the order of the halfedges incident to it. For each edge, it negates the associated curve. 
For each face, it inverts the order of the halfedges along its outer boundary. Similar to the 
overlay operation (see Section !^. 4. 21) . where the user can provide a set of ten functions, which 
are invoked when new vertices, edges, and faces of the resulting arrangement are created, 
while the overlay operation progresses, the user can provide a set of three functions that are 
invoked when a new vertex, halfedge, and face are created, while the reflection operation 
progresses. Extended data associated with these types, such as a primal vertex associated 
with an arrangement face as in the case of an arrangement representing a Gaussian map, 
can easily be updated with the provision of an appropriate function. 

The trivial observation that P © {—Q) = ~{{~P) © Q) leads to another optimization. 
Instead of reflecting all sub-parts in the set {P^ \ i = 1,2, . . . ,n, k = 1,2, . . . , rrii}, we reflect 
only the sub-parts in the set {P^ \i = 2, . . . ,n,k = 1,2,..., rrii}, and compute only the 
pairwise sub-parts Minkowski sums in the set {M^-^ 1 1 < i < j < n,k = 1,2, ... ,mj,i = 
1,2, ... , rrij}, their central projection, and the union of the appropriate projections to yield 
the set {Qij \ l<i<j<n}. Then, we apply the reflection operation described above 
on each member of this set, and obtain the full set of projections {Qij \ i = 1,2, . . . ,n, j = 
1,2, .. . ,n}. The Boolean flag associated with a face of an arrangement that represents a 
central projection is equal to the flag associated with its reflection. In other words, a face 
of Qij consists of directions that pierce Mij, if and only if, its reflection in Qji consists of 
directions that pierce Mjj. 

Phase 8 is purely topological. Thus, we do not expect the time consumption of this phase 
to dominate the time consumption of the entire process for any input. Nevertheless, it might 
be possible to reduce its contribution to the total time consumption through efficient testing 
for strong connectivity applied to all the DBGs [KMW98] . exploiting the similarity between 
DBGs associated with incident cells. Recall, that the set of arcs in a DBG associated with 
a vertex v is a subset of the set of arcs associated with an edge incident to v. Similarly, the 
set of arcs in a DBG associated with an edge e is a subset of the set of arcs associated with 
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Table 5.2: Time consumption (in seconds) of the execution of the eight phases applied to the Split Star 
assembly as input. Each one of the three rows refers to a different decomposition of the assembly. A 
— number of convex sub-parts per part. B — number of sub-part vertices per part. C — total number 
of convex sub-parts. D — total number of Minkowski sums. E — total number of arrangements of 
geodesic arcs embedded on the sphere constructed throughout the process. 



A 


B 


C 


D 


E 


1 


2 


3 


4 


5 


6 


7 


8 


3 


16 


18 


270 


607 


NA 


0.01 


0.04 


2.38 


0.41 


2.05 






5 


22 


30 


750 


1591 


NA 


0.01 


0.05 


5.03 


1.09 


7.07 


0.36 


0.01 


8 


32 


48 


1920 


3967 


NA 


0.01 


0.06 


11.12 


2.41 


27.99 







a face incident to e. The proposed technique reduces the cost from Oin^) per DBG to an 
amortized cost of 0{n^'^'^^), where n is the maximum number of arcs in any blocking graph. 

5.5 Experimental Results 

Our program can handle all inputs. However, we limit ourselves to a representative set of 
test cases, where we compare the impact of different decompositions on the process time 
consumption. The results listed in Table 15.21 were produced by experiments conducted on 
a Pentium PC clocked at 1.7 GHz. In all three test cases we use the Split Star assembly 
as input. Naturally, in all three cases identical projections are obtained as the intermediate 
results of Phase 6, hence the identical time consumption of the succeeding last two phases. 
Evidently, it is desired to decompose each part into as few as possible sub-parts with as small 
as possible number of features. However, an automatic decomposition operation may require 
large amount of resources to arrive at optimal or near optimal decompositions. Notice that 
Phases 4 and 6 dominate the time complexity. This is due to the large number of geometric 
predicates that must be evaluated during the execution of the overlay operation. 



Seriousness is the only- 
refuge of the shallow. 

Oscar Wilde 




Conclusion and Future Work 



In this thesis we show how a complete implementation of extendible arrangements with a rich 
set of operations enables a broad spectrum of robust applications that solve problems arising 
in domains such as motion planning, assembly planning, and solid modeling. For example, 
we describe how arrangements embedded on two-dimensional surfaces can be efficiently used 
to compute Minkowski sums of two polytopes in M^, which in turn, and in conjunction with 
several other operations based on such arrangements, can be used to partition an assembly 
with an infinite translation motion. The rest of this chapter is devoted to future prospects 
related to our research topics. 

6.1 Arrangements on Two-Dimensional Surfaces 

Constructing Minkowski Sums of polytopes in has been successfully attempted in the past. 
We introduce a robust, yet efficient method. Table 13751 shows that both our exact methods 
outperform the other exact methods. However, we believe that both of our methods, and in 
fact all Cgal based methods have great potential for further improvements through future 
optimizations applied to the infrastructure of Cgal, as Cgal is an evolving project. While 
the space consumption of the Cgm method is greater than the space consumption of the 
spherical Gaussian-map method, the table also reveals that the Cgm method is currently 
more efficient than its rival. We estimate that the gap will decrease, if not vanish, once all 
optimizations for the Arrangement_on_surf ace_2 data-structure that are still pending are 
implemented and enabled. 

We are constantly striving to improve the quality of our infrastructure, that is the 
Arrangement_on_surf ace_2 package. We have already identified few weak spots. Elimi- 
nating them will increase the genericity, extendibility, efficiency, and functionality of the 
package. We provide one example in each category. 
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6.1.1 Generic Observers 

There is a certain similarity between observers (see Section I2.4.4I) and visitors (see Sec- 
tion [23]T]), as typically each of their methods is triggered as a response to a certain event 
— a member of a pre-determined list of events. Technically, the main difference between 
them is that observers define a one-to-many mapping between objects, while visitors define 
a one-to-one mapping.^ Recall, for example, that a single arrangement may register many 
observers, but it is only natural to relate a single visitor to a specific algorithmic framework 
in order to realize a certain concrete algorithm. Consequently, arrangement observers are 
derived from a common base class, and their methods must be virtual. This is how mod- 
ules, which are closed for modification, are extended using object-oriented programming. 
However, composability of such modules is limited, since independently produced modules 
generally do not agree on common abstract interfaces from which supplied types must in- 
herit. In addition, when techniques from the object-oriented programming and the generic 
programming paradigms are mixed, they often clash. There are known methods to replace 
lists of objects, derived from a common base class, and linked during run time, with a list 
of syntactically unrelated objects concatenated during compile time (coded using a generic 
programming technique) [AleOll Chapter 3]. Nevertheless, we would like to simultaneously 
enjoy the benefits of both the object-oriented and the generic programming paradigms, that 
is, to enable the immediate production of composable modules that support dynamic poly- 
morphism. A very important research direction in our context is to explore these possibilities, 
perhaps pushing the limits of the C++ programming language along the way. 

6.1.2 Property Maps 

In many cases we need to associate values (called "properties") with the vertices, the 
halfedges, and the faces of the arrangement. In addition, it is often necessary to associate 
multiple properties with each vertex, edge, or face; this is what Boost literature refers to 
as multi-parameterization. BGL |SLL02] graph classes have template parameters for vertex 
and edge "properties" . A property specifies the parameterized type of the property and also 
assigns an identifying tag to the property. 

There are various ways to associate properties with arrangement cells. One option is 
to extend the geometric types of the kernel, as the kernel is fully adaptable and exten- 
sible |HHK"'"07] . However, this indiscriminating extension may lead to an undue space- 
consumption, as every geometric object is extended, regardless of its use.^ Another option, 
is to extend the vertex, halfedge, or face records of the DcEL; see Section [2.3.3[ This may also 
lead to excessive space-consumption, for example, when the data associated with a halfedge 
is in fact tied to the embedded geometric curve. In this case the data, or at least a handle 
to the data, must be stored twice in both twin halfedges. A third option is to extend the 
curve (or point) types defined by the geometry-traits class; see Section 12.3.11 But even this 
option leads to unjustified space-consumption, when only a limited number of arrangement 

^They also differ in essence. Wliile an observer typically implements a notifier, a visitor is usually a 
coherent part of an algorithm based on a fundamental and flexible framework. [GHJV95] 

^It also requires nontrivial knowledge about the kernel structure and the techniques to extend it. 
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features are associated with real data. In such cases it is advantageous to use external search 
structures that map individual arrangement features with their data. 

Designing a useful and convenient interface, while taking all considerations above into 
account, is a research topic on its own, which may further push the limits of good usage of 
the C++ programming language. 

6.1.3 Point Location for Surfaces 

Point location is one of the most fundamental operations applied to arrangements; see Sec- 
tion [2A51 The contest between the different point-location strategies for arrangement em- 
bedded in the plane was settled in favor of the "landmark" variants for many types of 
arrangements [HH08j . Unfortunately, at this point, these strategies cannot be applied on ar- 
rangements embedded on surfaces other than the plane. This is due to limitations that arise 
from the specific implementations of geometry traits that support these types of embedded 
surfaces, e.g., the geometry traits that handles geodesic arcs embedded on the sphere; see 
Section 12.61 The problem should be attacked from both its ends. That is, we can try to 
enhance the geometry traits implementations, and add all ingredients required by the con- 
cept ArrangementLandmarksTraits-2 as defined today, and at the same time, try to come 
up with alternative, perhaps similar, strategies that induce different requirements that are 
easier to satisfy by the geometry-traits classes. 

6.1.4 Geometry- Traits Models 

Continuous and steady effort is made to further extend the arsenal of geometry-traits models, 
or simply to improve the existing ones; see Section 12.5.21 Supporting arrangements induced 
by rich families of curves opens the door for numerous applications. 

The dominant bottleneck of all applications mentioned in this thesis is the application 
of the geometry operations implemented in the geometry-traits classes. Expediting their 
performance, while containing the growth of their memory footprint is always desired. For 
example, the arrangement package provides two traits classes that handle line segments. 
Both are parameterized by a geometric kernel; see Section [1.3.31 Segments defined by most 
Cgal kernels are represented only by their two endpoints. When a segment is split several 
times, the bit-length needed to represent the coordinates of its endpoints may grow exponen- 
tially (see |FWH04j for a discussion), which may significantly slow down the computation. 
Therefore, one of the two traits classes represents a segment by its supporting line in addi- 
tion to its two endpoints. When the traits class computes an intersection point of two line 
segments, it uses the coefficients of their supporting lines. When a segment is split at an 
intersection point, the underlying line of the two resulting sub-segments remains the same, 
and only its endpoints are updated. This traits class thus overcomes the undesired effect of 
cascading intersection-point representation, as described above, at the account of a larger 
memory footprint. A similar idea can be applied to the traits class that handles geodesic arcs 
embedded on the sphere. An implementation of a geometry-traits that handles such arcs 
that stores the projections of all the arrangement geometric features once calculated, and 
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retrieves them when subsequently needed, has great potential to reduce time consumption 
at the price of growth in space consumption. 

Another direction is to expand existing implementations to meet the requirements of 
the various concepts in the hierarchy described in Chapter [2J Consider, for example, the 
geometry traits class that handles geodesic arcs embedded on the sphere. Currently, it 
supports the basic operations required to construct and maintain arrangements induced by 
such arcs. As mentioned in the previous section, it might be possible to enhance it to enable 
the use of one or more of the variants of the "landmark" point-location strategies. Similarly, 
enabling Boolean set operations of spherical patches bounded by geodesic arcs embedded on 
the sphere requires the provision of few additional operations by the traits class. 



6.2 Three-Dimensional Arrangements 

Consider the following task: Given a set 5 = {^i, 5*2, ... , S'„} of two-dimensional surfaces 
in M^, construct the three-dimensional arrangement A{S) induced by S. Fulfilling this task 
in an efficient, complete, and robust manner has not been attempted yet, and is considered 
challenging. Implementing various strategies of point- location that operate on arrangements 
in and a plane-sweep and zone-construction frameworks for such a data structure is 
greatly desired, but extremely ambitious. In analogy to two-dimensional arrangements, a 
generic implementation of a plane-sweep framework can enable various operations, such as 
the overlay of spatial subdivisions and ordinary and regularized Boolean set operations of 
point sets bounded by general algebraic surfaces. These operations, in turn, can enable the 
implementation of a multitude of applications. 

Arrangements embedded on two-dimensional surfaces can be used as building blocks 
in the implementation of a data structure that represents a three-dimensional arrange- 
ment |Wein7j . We can consider each surface Sk separately, and construct the arrangement 
•A-k = -A-k^S) induced by intersection curves between Sk and S\ Sk embedded on S. The 
arrangements ^1,^2,- • • i-A-n can subsequently be connected together to properly represent 
the spatial subdivision A{S). 



6.3 Boolean Set-Operations 

In some sense and to some extent this thesis attempts to close gaps between theoretical 
results and practical needs. It is not accidental that great parts of the thesis are closely 
related to Cgal, as one of the goals of Cgal stated in the Introduction chapter of the thesis 
is to translate (theoretical) results into useful, reliable, and efficient programs for industrial 
and academic applications. Evidently, the Boolean set operations package of Cgal, which 
is based on the Arrangement_on_surf ace_2 package (see Section [2.7. ip is one of the most 
popular packages among Cgal packages in the commercial market. Naturally, we would like 
to continue improving this package. The problems addressed in the next two subsections 
were raised during the 3'"'^ Cgal User Workshop |dW08j [1]. 
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6.3.1 Fixing the Data 

Input data of Boolean set operations, namely, a set of one or more polygons, used in real- 
world applications is occasionally corrupted, as it originates from measuring devices that 
are susceptive to noise and physical disturbances. In some other cases, it contains many 
degeneracies, which either disable computations based on fixed-precision arithmetic, or slow 
down further computation using exact geometric computation. 



Invalid Data 




Figure 6.1: A relatively 
simple polygon that is 
not simple, given by its 
boundary {a, 6, c, d, b, e}. 



A polygon P is said to be simple (or Jordan), if the only points of 
the plane belonging to two edges of P are vertices of consecutive 
edges P P3|. Namely, no two edges intersect, except for every two 
consecutive edges, which share one endpoint. A simple polygon is 
topologically equivalent to a disk. A polygon P is said to be weakly 
simple, if the chain of the edges of P does not cross itself. A polygon 
P is said to be relatively simple, if it is weakly simple and the edges 
of P do not intersect in their relative interior. Observe, that a 
relatively simple polygon, the vertices of which appear only once in 
the boundary (the degree of each vertex is two), is simple. 

Input data for any Boolean set operation represents points set 
that may be bounded or unbounded, and may have holes. Items of 
such input take the form of general polygons or general polygons 
with holes with well-defined interiors and exteriors. A valid polygon 
must be weakly simple (but not necessarily simple) and its vertices 
must be ordered in counterclockwise direction around its interior. A 
valid polygon with holes that represents a bounded point set, has an 
outer boundary represented as a weakly simple (but not necessarily 
simple) polygon, the vertices of which are oriented in counterclock- 
wise order around its interior. In addition, the set may contain 
holes, where each hole is represented as a simple polygon, the ver- 
tices of which are oriented in clockwise order around the interior 
of the hole. Note that an unbounded polygon without holes spans 
the entire plane. Vertices of holes may coincide with vertices of the 
boundary. 

As mentioned above, real-world data is often corrupted. Naturally, passing invalid poly- 
gons (polygons with holes, respectively) as input to a Boolean set operation must be avoided. 
Apparently, automatically "fixing" corrupted data, that is, converting invalid input polygons 
or polygons with holes to valid ones, is not a simple task. Consider, for example, the self 
intersecting star depicted in Figure I^TST a). A point is considered inside the point set, if and 
only if the number of counterclockwise turns the oriented boundary makes around the point, 
also called the winding number, is greater than zero. It can be efficiently calculated using 
an Arrangement_2 data structure as follows. We extend each halfedge h with a Boolean flag 
that indicates whether the winding number increases or decreases when we cross h, that is. 




Figure 6.2: A poly- 
gon with a hole given 
by its outer boundary 

{a,b,c,d,e,f,g,h} and 
its hole {h, f, d, b}. 
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Figure 6.3: (a) A self crossing polygon given by {a,b,c,d,e}. (b) The Arrangement_2 data structure 
constructed from the polygon edges, (c) The Arrangement _2 data structure with updated face winding- 
numbers, (d) The Arrangenient_2 data structure with internal edges removed. 

when we move from the face incident to h to the face incident to the twin of h. Observe 
that twin halfedges always have opposite flag values. We extend each face with an integer 
that counts the winding number of every point in the face. Finally, we apply a Breath-First 
Search (BFS) on all the arrangement faces starting from the unbounded face and updating 
the face counters as we cross halfedges. The BGL, for example, can be employed for this 
task. Figure 16.31 illustrates the process. Once the winding number of every face is updated, 
we remove internal edges, as they are redundant, and convert the arrangement back into a 
valid polygon or polygon with holes. 

The problem becomes more complicated when the input polygons or polygons with holes 
violate several validity properties at the same time. Naturally, converting corrupted data 
into valid data consumes time. The challenge is to perform this task flawlessly and efficiently, 
while presenting a convenient interface to the user. 



Degenerate Data 

In computational geometry there are two main techniques to eliminate degeneracies and near- 
degeneracies. One is snap rounding [ GM95| IGGHT97t IHob99] and the other is controlled 
perturbation |HL04| lMO06j . Both techniques aim at processing geometric data, e.g., curves 
of the boundaries of general polygons, to yield new data that can be further robustly and 
more efficiently processed, perhaps using only limited precision. Traditionally, snap rounding 
has been applied to linear objects embedded in the plane |HP01l IHP02] . where it replaces 
sets of linear segments with sets of polylines. It can be extended though to other types of 
curves in the plane, such as Bezier curves |EKW07] . or even other curve types embedded 
on other surfaces that have a well defined grid (and perhaps other properties), such as 
geodesic arcs embedded on the sphere. Controlled perturbation has even a larger spectrum 
of applicable platforms. With respect to our polygons, applying any one of the two techniques 
above may result with (partially) overlapping curves, which belong to two different polygons, 
respectively. In this case, we need to merge the incremental winding contributions of the 
original curves mentioned above. 
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6.3.2 Improving the Efficiency 

The Boolean_set_operations_2 package provides efficient operations that compute the reg- 
ularized union or the regularized intersection of a set of input polygons. There is no restric- 
tion on the polygons in the set; naturally, they may intersect each other. The package also 
provides an efficient predicate that determines whether all polygons in a given set intersect. 

There are at least three different methods to compute the union 
of a set of polygons Pi, . . . Pm- We can do it incrementally as fol- 
lows. At each step we compute the union of Sk-i = [JlZi Pi with Pk 
and obtain 5*^. A second option is to use a divide- and-conquer ap- 
proach. First, we divide the set of polygons into two subsets. Then, 
we compute the union of each subset recursively, and obtain the 
partial results in Si and 5*2, respectively. Finally, we compute the 
union S*! U 5*2. A third option aggregately computes the union of all 
polygons. We construct an arrangement inserting the polygon edges 
at once, utilizing the sweep-line framework, (see Section 12.4.11) and 
extract the result from the arrangement. Similarly, it is also possi- Figure 6.4: The union 
ble to aggregately compute the intersection fli^i Pi '^^ ^ set of input of eight discs, 
polygons. 

The incremental method is more efficient for a small (constant) size of input polygons, 
and the aggregate method is more efficient for sparse polygons with a relatively small number 
of intersections. It is also possible to mix between the three methods, reaping the benefits 
of them all. We would like to figure out what are the exact conditions that should be used 
to determine when to use each method, or when to switch from one to another. 

6.3.3 Non Regularized Operations 

The Cgal package Nef_2 supports ordinary set-operations on point sets in |See07] . The 
point-set operands and results are rectilinear polygonal model. Such a point set can be 
defined by a finite set of open halfspaces, or obtained by set complement and set intersection 
operations. The package supports operations that consume and produce linear polygons 
defined by linear edges. The Boolean set operations package, on the other hand, supports 
only regularized set-operations, but the operations consume and produce general polygons. 
Recall, that a general polygon is a point set in that has a topology of a polygon, but 
its boundary edges map to arcs of curves, which are not necessarily linear. Extending the 
package to support not only regularized operations, but also ordinary ones, will make it 
useful for more applications; see, for example, Section 15.3.61 

6.3.4 Operating in 3-Space 

Boolean set operations are intuitive and therefore popular in many fields. For example, CSG 
is a representation model for solids based on Boolean set operations. Solids represented using 
CSG result from Boolean set operations applied to elementary solids called primitives, e.g., 
cubes, spheres, cones, and cylinders. A CSG solid is represented in a tree structure, where 
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the leaves represent primitives, and internal nodes represent Boolean operations. 

The Cgal package Nef_3 supports ordinary set-operations on point sets in |HK07j . 
Similar to the planar case, the package supports operations that consume and produce (lin- 
ear) polyhedra defined by (linear) halfspaces. Having the ability to construct and maintain 
arrangements in M^, will enable the development of a new package, that will support either 
regularized, or even non-regularized, robust Boolean set-operations that consume and pro- 
duce general (curved) polyhedra in M^, the boundaries of which are general surfaces. Many 
fundamental problems in solid modeling, motion planning, and other domains, can benefit 
from such a package. 

6.4 Collision Detection 

One possible progression of the collision detection algorithm and its implementation de- 
scribed in Section 13.41 is a complete integrated framework that answers proximity queries 
about the relative placement of polytopes that undergo rigid motions including rotation. 
The framework may use either the spherical or the cubical Gaussian-map to represent poly- 
topes. The interface of these two data structures should be consolidated to allow rapid 
interchanging. 

Some of the methods we foresee compute only those portions of the Minkowski sum that 
are absolutely necessary, making our approach even more competitive. Briefiy, instead of 
computing the Minkowski sum of P and — Q, we walk simultaneously on the two respective 
Cgm's, producing one feature of the Minkowski sum at each step of the walk. Such a strategy 
could be adapted to the case of rotation by rotating the trajectory of the walk, keeping the 
Cgm of —Q intact, instead of rotating the Cgm itself. 

6.5 Reflection Mapping and GIS 

We have developed a new data structure that can be used 
to construct and maintain cubical Gaussian-maps and com- 
pute Minkowski sums of pairs of polytopes represented by the 
new data structure; see Chapter [31 The name of the data 
structure is Cubical_gaussian_map_3, and we are considering 
introducing a package by the same name to a prospective re- 
lease of Cgal. The implementation is generic and can be 
used for other purposes, where six planar subdivisions embed- 
ded on a unit cube and stitched properly at the edges of the 
cube is useful, for example Cubical Environment Mapping; see, 
e.g., |FvDFH95[ Section 16.6]. 

In computer graphics, reflection mapping is an efficient 
method of simulating a complex mirroring surface by means 
of a precomputed texture image. The texture is used to store 
the image of the environment surrounding the rendered ob- 




Figure 6.5: Environment 
of the St. Peters Cathedral 
mapped on a teapot with a sil- 
ver material applied. Taken in 
RTHDRIBLE [22]. 
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ject. The surrounding environment can be represented, constructed, maintained, and stored 
in several ways; the most common ones are the Spherical Environment Mapping in which 
a single texture contains the image of the surrounding as reflected on a mirror ball, or the 
Cubical Environment Mapping in which the environment is unfolded onto the six faces of a 
cube and stored therefore as six square textures. 

Reflection Mapping can be categorized as some sort of a geographic information system 
(GIS). There are two broad methods used to store data in a GIS: Raster and Vector. In a GIS 
data is often related from different sources possibly of different storing types. Regardless of 
whether a single arrangement embedded on a sphere, or six arrangements embedded on the 
cube, are concerned, the connection to GIS is clear — both data-structures can accommodate 
geographic vector data in a natural way. 



6.6 Exact Complexity of Minkowski Sums 

Exact Maximum Complexity 



d 



d 



nil + 1712 + ... + nik 



Ei<i<i<fc(2"^i - 5){2mj - 5) + T.l<^<k^i + (2) 



The table to the right Dimension 
summarizes the known " 
exact bounds on the - 
maximum complexity of 
Minkowski sums of polytopes in terms of number of facets (((i— l)-faces) derived in Chapter|H 
The exact bounds are unknown for higher Dimensions as far as we know. 

It is known that the exact complexity (counting faces) of the Minkowski sum of two 
polytopes with m and n facets can be as low as m, when the two polytopes have the same 
number of facets m and parallel features, but it is unknown what is the minimum exact 
complexity of Minkowski sums of polytopes that have only a limited number of parallel 
features, or none at all. 
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Software Components, Libraries, and Packages 



A.l Visual Simulation 

We have developed a toolkit, called Sgal (Scene Graph Algorithm Library)/ that supports 
the construction and maintenance of directed acyclic graphs that represent scenes and models 
in M^. The toolkit includes, among the other, two interactive 3D applications. The first 
detects collisions and answers proximity queries for polytopes that undergo translation and 
rotation. The second enables users to visualize a scene in an interactive manner. It parses 
input files that describe the scene in a degenerate yet extended VRML format [27]. The 
format is degenerate, as not all VRML features are supported (yet). However, it has been 
significantly extended as described below. 

Both applications are linked with (i) Cgal, (ii) a library that provides the exact rational 
number-type, (iii) several Boost libraries, (iv) imagemagick [15] — a library that creates, 
edits, and composes bitmap images, and (v) internal libraries that construct and maintain 
3D scene-graphs, written in C++, and built on top of OpenGL. The internal code is divided 
into two libraries; SGAL — The main 3D scene-graph library and SCGAL — Extensions that 
depend on Cgal. 

We added several geometry nodes that represent polytopes using various sub-representations, 
such as Gaussian maps, a few other related nodes that handle coordinates, and many other 
miscellaneous nodes that provide services, such as antialiasing and snapshoting. The de- 
scriptions of some of the geometry nodes follows. 

ArrangementOnSphere This node represents models as arrangements of geodesic arcs em- 
bedded on the sphere. 

ExactPolyhedron This node represents polyhedra as meshes using the Cgal Polyhedron_3 



We plan to offer Sgal with an open-source license in the future, making it available to the public. 
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data structure. 

SphericalGaussianMap This node represents polytopes as spherical Gaussian maps using 
the Arrangement_on_surf ace_2 data structure instantiated as an arrangement embed- 
ded on the sphere. 

CubicalGaussianMap This node represents polytopes as cubical Gaussian maps using the 
Cubical_gaussian_map_3 data structure. 

The implementation relies on inputing exact coordinates. To this end, the format was 
further extended with a node called ExactCoordinate that represents exact coordinates, 
and enables the provision of exact rational coordinates as input. 

The entire partitioning process described in Chapter [5] is realized within Sgal. The 
library contains all the necessary ingredients required to represent and visualize the input 
and the output, and to simulate the process. In particular, it has been extended with two 
geometry node types: the Assembly node type represents assemblies or subassemblies, and 
the AssemblyPart node type represents parts of assemblies. Notice, that each node ob- 
ject of the three types AssemblyPart, SphericalGaussianMap, and ArrangementOnSphere 
internally maintains the Cgal data structure that represents an arrangement of geodesic 
arcs embedded on the sphere |BFH+07] . an instance of the Arrangement_on_surf ace_2 class 
template. 



A. 2 Software Availability 

Compiling and executing the programs require the latest internal release of Cgal (post 
version 3.31) and the components listed above including an internal package of Cgal that 
supports the Cubical_gaussian_map_3 data structure. The source code is available upon 
request.^ Precompiled executables compiled either with g++ 4.2.3 on Linux or with VC 8 on 



Windows, data sets, and documentation can be downloaded from |http: //www, cs . tau. ac . 



il/~efif/CD/3d, 



^Send email to ef if ogeKSgmail .com . 
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